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I.  INTRODUCTION 


The  existence  of  flaws  in  solid  propellant  grains  is  a  problem  that 
has  plagued  the  solid  rocket  motor  industry  for  a  considerable  time/1’2,3^ 
Flaws  may  exist  in  the  grain  due  to  unsteady  processes  during  casting, 
difficulties  in  curing,  problems  in  humidity  and  temperature  control,  etc. 
The  surveillance  programs  that  have  been  under  way  for  several  years 
indicate  that  flaws  may  form  due  to  "aging. 

The  Air  Force  Rocket  Propulsion  Laboratory  at  its  cwn  laboratory  and 
through  its  contractors  has  carried  out  numerous  investigations 
dealing  with  the  problem  of  initiation  of  burning  within  flaws  ii  solid 
propellant  grains.  Similar  studies  have  been  carried  out  in  Japan ^  and 
the  United  Kingdom. ^  The  tests  conducted  by  AFRPL  and  its  contractors 
primarily  dealt  with  ignition  times  and  rates  of  propagation  of  the  flames 
into  "manufactured  cracks"  although  tests  were  conducted  at  Tniokol's  Utah 
Division  in  which  the  fracture  was  obtained  by  tearing  the  propellant.  The 
tests  conducted  in  Japan  were  also  conducted  with  manufactured  cracks.  The 
United  Kingdom  tests  conducted  at  the  Explosive  Research  Development  Estab¬ 
lishment  were  tests  in  small  scale  rocket  motors  with  manufactured  cracks 
and  debonds. 

Tests  have  been  conducted  at  Thiokol,  Utah  Division^  and  at  the 
Naval  Weapons  Center,  China  Lake ^ 10 ^  dealing  with  debonds.  These  tests 
were  all  conducted  utilizing  a  'half -motor;"  Hie  "half-motors"  were  half 
of  a  solid  propellant  motor  and  a  block  of  plexiglas  joined  together  to 
allow  visual  observation  of  the  burning  processes.  Similar  unpublished 
tests  have  been  conducted  by  H.  R.  Jacobs  at  the  University  of  Utah. 

All  of  the  above  mentioned  tests  dealing  with  burning  in  cracks  and 
debonds  indicated  overwhelming  statistical  evidence  that  burning  would 
propagate  into  cracks  and  debonds  and  that  the  propagation  rates  increased 
with  the  external  or  chanter  pressure.  In  some  of  the  half-motor  tests 
catastrophic  failure  occurred. 


The  present  study,  which  was  the  result  of  an  unsolicited  proposal, 
had  as  its  primary  objective  the  determination  of  whether  the  pressures 
which  develop  in  burning  flaws,  cracks  and  debonds  are  sufficient  to 
cause  the  defect  to  propagate  due  to  mechanical  failure  of  the  propellant. 

The  University  of  Utah  in  its  investigation  was  directed  to  assume 
that  the  defect  surface  was  ignited.  With  this  assumption,  it  was 
further  assumed  that  the  size,  shape  and  nature  of  the  initial  defect 
were  known. 

In  order  to  carry  out  the  proposed  study  a  qualitative  investigation 
of  the  governing  parameters  was  first  conducted.  The  critical  parameters 
are,  in  addition  to  propellant  characteristics,  main  chamber  pressure  and 
defect  geometry,  the  main  chanter  gas  velocity,  the  gas  velocity  within 
the  defect,  the  burning  rate  of  the  wall,  and  the  fracture  propagation 
velocity. 

If  the  burning  rate  of  the  material  is  greater  than  the  crack  propa¬ 
gation  velocity,  the  defect  tip  will  be  "burned  out"  before  it  can  be 
propagated  due  to  mechanical  failure.  If  the  surface  burning  of  the  defect 
area  is  large  compared  to  its  cross-sectional  area  at  its  exit  plane,  then 
the  velocity  of  the  burned  gaseous  effluent  will  in  general  be  larger 
than  that  external  to  the  defect  and  there  will  be  no  net  flow  imo  the 
defect  but  only  a  flow-out.  The  most  important  flow  is  that  of  the 
effluent  gases,  for  it  is  this  velocity  whic1'  determines  the  pressure 
field  within  the  flaw.  For  most  solid  propellants  the  combustion  rate  is 
pressure  dependent.  The  stress  field  around  the  defect  is  also  dependent 
upon  the  pressure,  and  hence,  the  presence  of  "cracking"  and  its  propagation 
speed. 

The  pressure  distribution  in  a  burning  defect  will  vaiy  with  time 

because  the  defect  geometry  is  being  continually  changed  by  burning  at 

the  walls.  In  addition  cracking,  if  it  occurs,  will  alter  the  geometry 

and  thus  change  thepressure  distribution.  Despite  these  facts,  a  true 

transient  analysis  is  not  necessary  if  the  geometry  change  is  slow 

compared  to  the  velocity  of  the  combustion  products  within  the  flaw.  In 

general,  the  propellant  burning  rate  is  on  the  order  of  10  *  fps  while 

2 

that  of  the  gas  velocity  is  on  the  order  of  10  fpsj  therefore,  in  regard 
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to  the  burning  wall,  a  quasi-steady  analysis  is  justifiable.  The  quasi¬ 
steady  analysis  as  related  to  crack  propagation  speed  is  reasonable  if  the 
crack  propagation  velocity  is  either  much  slower  than  the  gas  velocity  or 
if  the  defect  remains  stable  for  a  length  of  time  and  then  propagates  to  a 
new  stable  geometry  which  is  instantaneously  ignited.  With  these  facts  in 
mind,  the  current  study  has  limited  the  gas  phase  study  tc  a  quasi -steady 
analysis  to  determine  the  instantaneous  pressure  distribution  in  a  burning 
defect  in  a  solid  propellant  grain. 

Primarily,  the  analyses  to  determine  pressure  distributions  have 
been  limited  to  one  dimensional  analyses.  These  analyses  are  applicable 
primarily  to  flaws  with  large  burning  surface  area  compared  to  the  exit 
flow  area. 

A  two-dimensional  analysis  was  carried  out  to  determine  the  effect  of 
two  dimensions  on  triangular  cracks  with  small  angles  of  divergence.  The 
results  of  this  analysis  which  have  been  included  in  Section  II  justifies 
the  one-dimensional  assumption. 

A  parametric  study  was  conducted  to  determine  the  effects  of  geometry 
and  fuel  characteristics  on  pressures  within  cracks  and  debonds  (burning  on 
one  side  of  the  crack  only)  and  is  reported  in  Section  III.  Also  included 
in  this  section  is  a  study  indicating  typical  times  required  for  the 
pressure  at  the  crack  tip  to  decrease  by  a  factor  of  two  and  indications  of 
crack  geometry  changes  with  time. 

Section  IV  of  the  report  sumnarizes  the  work  for  cracks  in  pressurized 
cylinders  conducted  to  determine  criticality  conditions  by  the  use  of  time- 
dependent  stress  intensity  factors.  Methods  of  determining  times  to  failure 
for  a  pressurized  crack  and  means  of  determining  initial  crack  propagation 
velocities  through  the  use  of  a  thermodynamic  power  balance  are  also  included. 

Section  V  compares  crack  propagation  velocity  estimates  to  the  propellant 
burning  rate  for  TPH-lOil  and  estimates  of  times  to  failure,  while  Section 
VI  summarizes  the  conclusions  of  this  program. 


II.  DETERMINATION  OF  PRESSURE  DISTRIBUTION 
IN  BURNING  FLAWS,  CRACKS  AND  DEBONDS 


ONE-  PI  MENS  TONAL  ANALYSTS  TNinjlBING  FRICTION 

The  relations  governing  the  state  of  the  flow  are  independent  of 
geometry  :iikI  are  stated  below: 


liquation  of  State 

P  =  pRT 

U) 

Energy 

h*u2/2  *  hQ 

(2) 

Heat  Capacity 

h-hr  =  CpT 

(3) 

In  the  energy  equation  it  is  assumed  that  the  mass  addition  is  accomplished 
at  constant  enthalpy  and  negligible  kinetic  energy.  This  a  reasonable 
assumption  for  a  burning  surface  if  the  combustion  zone  is  assumed  infinitely 
thin  and  at  the  surface.  It  is  also  assumed  that  the  flow  is  adiabatic. 

The  latter  is  a  good  assunption  for  solid  propellant  grains  which  release 
considerable  energy  and  are  notably  poor  conductors .  The  specific  heat  and 
viscosity  are  assumed  constant,  implying  that  no  chemical  reaction  takes 
place  after  the  gas  is  added  to  the  stream. 

Ihe  one- dimensional  momentum  equation  is 

7  p  af  *  axW  *  Twb  -w  -  Twn  ~cGr  u)  (4) 

The  change  of  y  with  x,  dy/dx,  is  assumed  small  to  insure  one-dimensionality 
and  since  normal  burning  is  assumed  the  x  component  of  momentum  of  the  small 
mass,  dm,  added  over  the  length  dx,  is  negligibly  small.  The  variation 
of  flow  area.  A,  and  differential  surface  area,  dS ,  with  x,  is  of  the  form 


y(x)x£ 

(5a) 

x*dx 

(5b) 

whore  the  exponent  £  may  assume  the  values  0  or  J.O  depending  on  how  the 
channel  depth  varies  with  x.  To  elucidate  this  point  consider  the  conical 
shaped  duct  illustrated  in  Figure  1. 
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Figure  1,  Conical  duct 

The  area  at  any  r#  is  given  by- 

A  *  2*yr  (6a) 

where 

r  =  ro  ’  x  (6b) 

and;  therefore,  this  geometry  corresponds  to  £  =  1. 

Three  types  of  geometries,  which  appear  frequently  in  application  and 
to  which  Equation  4  applies  are  illustrated  in  Figure  2. 

The  geometry  labeled  "crack"  represents  a  one -dimensional  model  of  a 
deep  flaw  within  the  propellant  grain  itself.  The  debond  geometry  is  the 
one-dimensional  model  of  a  separation  between  the  propellant  grain  and  the 
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cylindrical  portion  of  the  casing,  while  the  circular  deb and  is  the  model 
of  a  similar  separation  between  the  grain  and  the  hemispherical  bulkhead  of 
the  casing. 

Observation  of  Figure  2  yields  the  following  specific  forms  for  A  And  dS 

Crack  A  =  y 

dsn  .  o 

dS^  *  2dx 

Deband  A  =  y  (7) 

*  d^  «  dSn  -  dx 

Circular  Debard  A  -  2iry(ro  -  x) 

-  dS^  =  2ir(ro  -xjdx 

Ihe  fact  that  surface  area,  d£,  may  be  directly  related  to  distance  along  the 
flaw,  dx,  follows  from  the  assumption  that  dy/dx  is  small. 

Ihe  general  momentum  equation  contains  wall  shear  stress  terms  which 
must  be  related  to  more  readily  measurable  quantities.  Two  tvpes  of  wall 
shear  stress  are  encountered;  shear  at  a  burning  wall,  7^,  and  shear  at 
a  non-burning  wall,  t^.  Both  types  of  wall  shear  stress  are  assumed  to 
be  approximated  by  a  relation  of  the  fore 

Tw  =  Ifpu2  C8) 

The  friction  factor,  f,  fur  the  non-burning  wall  is  solely  a  function  of 
Reynold’s  number.  The  friction  factor,  f^,  for  the  burning  wall  is  dependent 
both  on  Reynold's  number  and  the  mass  addition  rate  since  this  mass  addition 
tends  to  "blow”  the  boundary  layer  away  from  the  wall  and  thus  reduce  wall 
shear  stress.  Reference  12  has  reported  a  correlation  between  friction  factor 
for  external  flow  with  pressure  gradient  and  non-dimensional  blowing  rate, 

Br,  where 


(9) 


B 

r 


and  Re^  is  the  Reynold’s  number  based  upon  the  distance,  x,  from  the 
leading  edge.  In  the  problem  at  hand, which  more  closely  resembles 
pipe  flow  than  boundary  layer  flow,  it  is  more  reasonable  to  base  the 
Reynold's  number  on  the  local  hydraulic  diameter. 

For  a  parallel  slot  the  characteristic  length  used  is  the  crack 
(or  debond)  height,  y. 


R 


e 


=  zyup 

v 


(10) 


Since  the  problem  under  consideration  is  one-dimensional  P,.r  =  p  =  p 

w  S 

and  the  normal  velocity  at  the  wall,  V  ,  can  be  evaluated  by  considering  a 
differential  wall  length  dx,  in  which  case  one  has 

.  dm  «  p  Vw  dx  (11} 


Utilizing  the  above  relation  in  Equation  9  one  obtains 


Br  "  ~  * 

oU 


dm 

dx 


C12) 


The  correlation  of  Reference  12  may  be  approximated  to  obtain  a 
friction  factor  with  mass  addition  of  the  foim 


fb  =  Hf  (13a) 

where 

H  '  7  [74x7  *  pS7T]  l13b> 

Equations  8,  13a,  and  13b  determine  the  wall  shear  stress  on  both  burning 
and  non-buming  walls  once  the  functional  relationship  between  the  friction 
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factor,  f,  and  the  Reynolds  nimber,  R^,  is  obtained 
is  of  two  basic  types  depending  on  whether  the  flow 
For  laminar  flow  the  relation  is 


This  functional  relation 
is  laminar  or  turbulent. 


f  = 


96 

IT 


(14a) 


and  for  turbulent  flow 


0.316 


(14b) 


Substituting  Equations  7,  S,  and  13a  into  the  general  momentum 
equation,  one  obtains  the  specific  momentum  equations  for  the  geometries  of 
Figure  2  as 


Crack 


y§ -  ^  -  k(mu)  J 


(15a) 


/ 


P?bond 


dP  n  2 

xaj  - 


1  ♦  H 
2 


d  (  m  u) 
3x 


(15b) 


Circular 

Debond 


2*r[4  -  fei*)  - 


fpu 


2  1  +  H- 


d  (  k  u) 

33T 


(15c) 


where  f  and  H  are  determined  by  Equations  12,  13b,  14a  and  14b.  Notice 
that  the  momentum  equation  for  the  circular  debond  Equation  15c  reduces  to 
the  same  equation  as  that  for  the  debond  Equation  15b  if  one  evaluates  the 
limit  as,  r  -*•  <*>. 

One  other  equation  determines  the  motion  of  the  flow;  this  is  the 
continuity  equation 

m  =  puA  (16) 


For  any  given  geometry  there  are  two  equations  of  motion  and  three 
equations  of  energy  and  state.  These  equations  are  expressed  in  terms  of 
seven  variables  h,  T,  y,  u,  P,  p  and  m.  Since  there  are  five  equations 
and  seven  variables,  it  is  possible  to  express  any  one  parameter  as  a  function 
of  any  other  two. 

In  general,  experimental  data  are  available  to  describe  the  burning  rate 
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of  the  material.  Since  all  mass  flow  in  the  flaw  is  due  to  the  burning 
rate  which,  in  general,  is  dependent  on  the  pressure,  and  since  the  original 
geometry  of  the  flaw  is  a  known  function  of  x,  it  is  possible  to  express 
P  as  a  function  of  m  and  y.  The  simultaneous  solution  of  the  governing 
equations  is  initiated  by  using  the  perfect  gas  relation  Equation  IV  the 
property  relation  Equation  3  and  the  continuity  relation  Equation  16  to 
eliminate  p,  u,  and  h  from  the  momentum  and  energy  equation  yielding 

jfy'  (V  *  V  -  -  5  T"ZR<sn’  *  “V  =  s4t"2r>  (17a) 

(y  -  l)(Tn2R)2  +  (2P2Y)(Tn2R)  -  2TQ  F2Y  n2!*  =  0  (17b) 

(the  prime  denotes  differentiation  with  respect  to  x) 
where 

F  =  P/Pr  and  n  =  m  /PjA. 

Pr  is  a  referenced  pressure  and  may  be  picked  at  any  desired  value. 

Equation  (17b)  is  new  solved  for  Tn2R  by  the  quadratic  formula  to  obtain 

tA  -  *.  gl^ 

Y  -1  (18) 

where  G  =  2(y  -  1)TqR/y  is  a  constant  parameter  of  the  propellant.  Price 
has  shown  that  of  the  two  alternate  signs  only  the  plus  sign  has  physical 
meaning  since  the  minus  sign  corresponds  to  removal  of  mass  in  an  unmixing 
process  which  leads  to  a  decrease  in  entropy.  He  further  defined  the  radical 
as  a  new  variable 

c  =  (F2  +  gY)*5 


Equation  18  is  now  combined  with  Equation  17a  and  rearranged  to  obtain 


dF 

dx 


[a'  -  f  (V 


+  U 


n  [!  <Sn 


AG2nn'  | 


(19) 


+ 
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To  proceed  further  a  mass  addition  mechanism  must  be  postulated. 

If  errosive  burning  is  neglected  the  combustion  rate,  rfa,  of  a  solid 
propellant  way  be  assumed  to  follow  the  empirical  relation 

rfa  =  CP11  (20) 

where  rfa  xs  the  recession  rate  of  the  burning  wall  and  C  and  n  are 
constants  to  be  determined  empirically  for  a  particular  fuel.  The 
local  differential  mass  flux  is  then  drti  =  Dr^,  where  D  is  the  density 
of  the  solid.  Hie  local  mass  flux  may  new  be  integrated  to  obtain 

A 

*  -J0  "Vft,  ♦  %  (21) 

The  flip  indicates  the  presence  of  burning  on  the  head  Kid.  The  value  of  ift 
is  0 


A  = 


,n. 


(22) 


Combining  Equations  21  and  22  with  the  definition  of  n,  one  obtains 


~r 


do  ^ 


n0  ■  M’o'1 


(23a) 

(23b) 


Determination  of  the  pressure  distribution,  or  equivalently  P,  as  a 
function  of  x  amounts  to  a  simultaneous  solution  of  the  differential 
equation  for  F  (Equation  19)  and  the  integral  equation  for  n,  Equation  23a/ 
For  the  three  geometries  of  interest  the  equations  to  be  solved  are 
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In  order  to  solve  the  governing  equations,  Equations  25,  26,  or  27,  it 
is  necessary  to  establish  the  boundary  condition.  If  the  flow  is  subsonic 
at  the  exit  plane  of  the  flaw,  then  the  pressure  must  equal  the  main 
chamber  pressure,  P^.  If  the  flow  chokes  at  the  exit  plane,  the 
boundary  condition  is  just  Mach  number  equal  to  one.  Although  it  is  pos¬ 
sible  for  flow  in  a  diverging  duct  with  mass  addition  to  reach  supersonic 
velocities,  it  is  not  a  normal  situation.  To  ascertain  whether  the 
exit  Mach  number  can  become  greater  than  one  it  is  necessary  to 
evaluate  the  derivative  of  the  Mach  number  as  it  approaches  one. 

Since  the  Mach  number  is  of  such  prime  importance,  it  should  be 
introduced  into  the  problem  formulation.  Because  a  perfect  gas  equation 
of  state  has  been  assumed  as  well  as  constant  specific  heats,  the  standard 
relationships  developed  in  Reference  14  are  applicable  for  Mach  numbers. 

'Hie  term  m/P^A  which  was  given  the  symbol  n  can  be  written  in  terms 
or  temperature  and  velocity 


n 


(28) 


2u 

W 
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The  temperature  may  be  elininated  in  Equation  28  utilizing  the  adiabatic 
expansion  relationship 


(29) 


and  the  velocity  u  may  be  eliminated  using  the  definition  of  Mach  number 
and  the  perfect  gas  formulation  for  the  speed  of  sound. 


After  rearranging  the  terms  there  is  obtained 

t2  .  n2  2,h 


m2  ,  ctr  ♦  cvr  -  f) 

F(y  -  1) 

Utilizing  the  previously  introduced  notation 


(30) 


?  -  tf  +  g2!,2]*5  , 


Equation  30  simplifies  to 
M2  - 

F(t  -  1) 


(31) 


In  reality,  only  the  plus  sign  in  Equation  31  has  significance  as  has 
been  pointed  out  by  Price  . 

The  solution  of  the  simultaneous  governing  equations  (Equations  25, 
26,  or  27)  may  new  be  determined  with  the  restrictions  placed  on  the 
boundary  exit  plane  of  either 


(a) 

Subsonic  flow 

P  =  P, 
e  ch 

(b) 

Choked  flow 

o 

• 

»-H 

II 

(c) 

Supersonic 

.  ..  dM2 

when  lim 

is  positive 

M  -*■  1 


The  condition  (c)  may  be  evaluated  utilizing  a  general  relationship 
developed  by  Shapiro. ^  Application  of  Shapiro’s  general  formulation 
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where 


M*(l +  1-^-1  M2  £ 


2  g  (s;  ♦  h^)  ♦  3L ;  **?  ^ 


From  Equation  33  it  is  clear  that  Aether  the  local  Mach  number 
increases  or  decreases  depends  on  whether  the  local  Mach  number  is 
greater  or  less  than  unity  and  also  whether  the  G(x)  is  positive  or 
negative.  Since  in  the  case  under  consideration  the  Mach  number  increases 
from  practically  zero  at  the  crack  tip,  the  sign  of  G  determines  whether 
locai  Mach  number  increases.  The  value  of  G  is  primarily  dependent  on 
the  area  and  mass  addition  functions  and,  therefore,  the  sign  of  G  is 
determined  by  the  crack  geometry  and  fuel-burning  properties.  For  typical 
fuels  and  geometries,  G  is  always  positive.  The  Mach  number,  therefore, 
increases  toward  unity  and  unless  G  happens  to  decrease  to  zero  precisely 
as  M  *  1  (which  in  general  is  not  true) ,  the  flew  will  choke  at  this 
point.  When  choking  occurs,  there  is  a  transient  period  of  readjustment 
in  which  the  tip  pressure,  Pq,  is  increased  by  an  amount  such  that 
Mach  1  occurs  at  the  exit  plane,  the  pressure  there  is  greater  than  the 
external  pressure  and  expansion  waves  form  at  the  flaw  exit  and  extend 
into  the  rocket  chamber. 

The  complicated  nature  of  the  governing  equations  together  with  the 
pressure  dependence  of  the  empirical  burning  laws  for  solid  propellants, 
in  general,  allow  for  only  numerical  solutions  of  the  stated  problem. 
Appendix  A  describes  the  numerical  procedures  which  were  used  to 
generate  solutions  to  the  problem  of  mass  addition  with  friction  in 
variable  area  ducts. 
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APPROXIMATE  ANALYTICS  SOLUTION  WITHOUT  FRICTION  FOR  SPECIAL  MASS  FLUX 
DISTRIBUTIONS 


An  analytical  solution  may  be  obtained  for  the  pressure  distribution 
i:.  a  crack  or  debond  if  one  neglects  the  effects  of  wall  shear  and  the 
mass  flux,  n,  varies  as  where  k  and  q  are  constants  depending  upon 
the  geometry  and  fuel  characteristics. 

Equation  19  with  the  friction  factor,  f,  equal  to  zero  reduces  to 

F'*  (t  -\p)  (lt  '  n  y"  *  ^f")  (W) 

for  cases  whe re  the  geometry  is  such  that  i  -  0. 

The  form  of  Equation  34  is  such  that  the  relationship  between  Gn} 

7  and  c  corresponds  to  the  sides  of  the  triangle  shown  in  Figure  3. 


Figure  3.  Relationship  of  Gn,  F  and  t. 

As  can  be  seen  from  Figure  3,  t,  F  and  P*  can  be  expressed  in  terms 
Gn  and  an  angle  4. 


F' 


c  =  -SSL 

COS<f> 

F  =  Gn  tan 


Gn 


ZoFJ*'  +  tan*  G* 


(35) 

(36) 

(37) 


Substituting  Equations  35,  36  and  37  into  Equation  34  and  simplifying, 
there  is  obtained 
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*'  -  T  ■-■‘yliM  (rU  '  si“*)  £  *  h  ’  sin  ♦>  f)  (38) 

If  y'/y  '  n'/n,  from  which  it  follows  y'/y  ^  m'/m,  and  if  y  and  m  are 
known  functions  of  x,  the  possibility  exists  for  separation  of  variables. 
A  particular  form  for  which  this  is  true  Is 


m  '  x  (39) 

y ' *q  (40) 

In  reality,  m  is  given  by  the  integral  equation.  Equation  21.  How¬ 
ever,  under  particular  situations  m  can  be  approximately  fitted  by  a  power 
function  such  as  Equation  39. 

Substituting  Equations  39  and  40  into  Equation  38  yields 

♦'  =  xH^TsinO  '  [Yk  '  CqCY  -  1)  +  k)sin*] 

Upon  separation  of  variables  there  is  obtained 

dx  _  (1  -  y  sinA)  d* 
x  B(k,  q,  *j  cos* 


(41) 


where 


B(k,  q,  4>)  =  rk  -  (q(y  -  1)  +  k)  sin* 


(42) 


(43) 


Equation  42  falls  into  one  of  three  distinct  types,  depending  on  the  re¬ 
lationship  of  k  to  q.  The  three  types  are  as  follows: 

1.  k  =  q,  therefore,  B(k,  q,  *)  =  Yk(l  .  sin^ 


^  . 


k  =  Ty^TT’  therefore,  B(k,  q,  *)  =  *(l  +  sin^ 


3.  k  ^  m  and  k  +  -q  jl 
n  (y  +  IJ 

If  the  relation  between  k  and  q  is  of  Type  1  or  2,  then  upon  integration 
of  Equation  42  one  obtains 


x 


ta4*$)  Yq  [;  7n  (1  Vs-inTr] 


(44) 
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where  +sin$  corresponds  to  k  *  -q  and  -  sin$  to  k  =  q.  If  the 

relation  between  k  and  q  is  of  Type  3,  then  Equation  43  integrates  to 


lk(>l)  «q(}-l]]Tq=Jc7 


(45) 


In  order  to  determine  what  values  of  <|iare  relevant  in  Equation  44  and  45, 
one  must  consider  the  Mach  number.  Applying  the  transforms  of  Equations 
35,  36  and  37  to  the  Mach  number  relation,  Equation  31,  the  following  is 
obtained. 

»2  _  1  -  sin « 

M  sin*  I7T  (46) 


From  Equation  46  it  can  be  seen  that  values  of  $  less  than  sin'1  1/y 
correspond  to  M  >  1.0  and  values  of  +  greater  than  sin"1  1/y  to  M  <  1.0. 
Since  the  governing  equation  was  derived  assuming  M  <  1.0,  it  is  apparent 
that  $  must  be  greater  than  or  equal  to  sin-1  1/y 

Considering  Figure  3,  it  can  be  seen  that  since  both  Gn  and  F  are 
positive  quantities,  $  must  be  less  than  or  equal  to  n/2.  Now  again 
considering  Equation  45,  notice  thft  the  tern  on  the  far  right  may  take 
on  imaginary  values  if  B(k,  q,  <|>)  becomes  negative.  The  possibility  of 
B(k,  q,  $)  taking  on  negative  values  depends  upon  the  relationship  of 
k  to  q  and  upon  the  values  idiich  4>  is  allowed  to  assume.  Setting  <f>  =  n/ 2 
and  B(k,  q,  <j>)  *  0  and  solving  for  k  in  terms  of  q,  one  obtains  k  =  q. 
Therefore,  if  k  is  greater  than  q,  the  term  B(k,  q,  <|>)  does  not  become 
negative  and  the  upper  limit  on  4  is  n/2.  If  k  is  less  than  q,  the 
term  B(k,  q,  <*>)  may  become  negative  and  in  this  case  the  upper  limit  on 
<j>  is  that  value  of  $  at  which  B(k,  q,  $)  =  0.  As  was  pointed  out  before, 
the  smallest  value  which  <j>  may  assume  is  sin-1  1/y  at  which  point  M  =  1.0; 
therefore,  setting  4>  -  sin_1l/yand  B  (k,  q,  <f>)  =  0  and  solving  for  k, 
one  obtains  k  =  q/(y  +  1).  This  last  result  indicates  that  k  must  be 
greater  than  or  equal  to  q/(y  +  1)  for  the  equation  to  yield  a  realistic 
flow.  It  is  apparent  that  the  +sin#  choice  in  Equation  44  has  ro  physical 
meaning  since  it  corresponds  to  a  flow  which  is  choked  at  i  =  9,  All  of 
the  above  results  are  summarized  in  Figure  4. 
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k  >  q 

sin”1  i  <  £  <  ir/2 

Y  ~  — 

k  <  q 

sin  -<  4  <  sm  5-+-q^"iy 

V  '  ^  . 

flow  choked 

-  -  Y  +  1 

Figure  4.  Dependency  of  relevant  <j>  range  on  the  relationship  of  k  to  q. 


For  the  relevant  range  of  <f>  enumerated  in  Figure  4,  one  may  consider 

the  relationship  between  <f>  and  x.  Notice  in  Equation  45  that  if  k  >  q, 

then  [k( y +  1)  +  q(Y  -  l)](q  -  k)  <  0  and,  therefore,  as  <t>  -*•  ir/2,  x  -*■  0. 

When  k<q,  x  =  0at<j>  =  sin  ^  ijc/(k  +  q(Y  -  1)).  If  k  =  q.  Equation  44 

describes  the  relationship  of  x  to  <f>.  Letting  $  -*•  ir/2  in  Equation  44, 
note  that  tan(ir/4  +  4>/2)  2/(1  +  cosif>  -  ri:4).  Setting  1  -  sin<|>  =  w  and 

letting  w  0  is  equivalent  to  letting  <p  ir/2;  therefore,  making  the  above 
substitution  in  Equation  44,  one  has  as  4>  v/2 


x  -*■  const 


1  +  Y 

2y<T 


1  -  Y 
2yqw 


But  the  above  limit  is  of  the  form 


x 


S. 

w 

const  77 
w 


1  +  Y 
2Yq 


where 

a  -  1  -  I 
g  "  1  +  y 


Expanding  eg//w  into  a  series  and  noting  that  g  <  0,  one  obtains  as  <j>  -*■  tt/2, 
x  ->  0  in  Equation  44. 
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From  the  definition  of  n  and  the  assumed  fora  of  y  and  ro  it  is 
apparent  that 


n  ^  ^ 


Taking  both  sides  of  Equation  45  to  the  k 
that  if  l  f  q 


q  power,  one  can  determine 


k(Y+i)|  +  i)  +  q(r 


If  k  =  q,  however,  n  is  a  constant. 

■me  parameter  of  greatest  concern  is  F  and  the  region  of  greatest 
concern  rs  the  tip  region.  Mien  ]c  >  q,  it  has  been  shown  that  x  »  0 
corresponds  to  ♦  ->  „/2.  action  43  gives  „  as  a  function  of  *  and 
Equation  36  gives  F  as  a  function  of  ,  and  ♦.  Therefore,  substituting 
Equation  48  into  Equation  36  gives  F  as  a  function  of  »  alone 


P  '  tan$ 


(tan  (j  *  ^ 


q  -  k(y  +  1)1 


(y  +  1J  +  q(y 


Observing  Equation  49  as  ^  „/2,  one  may  note  the  following 

q,  ♦)  -*•  (y-  l)(k  -  q) 

therefore,  as  $  tt/2 


F  •+  Constant 

cos$ 


L _ A _ p/_i_\ 

,  \1  +  COS <t>  -  sin<j>/  \cos <pj 


q-k(y+l)  (MY+TFqOr 


Simplifying  and  rearranging,  one  has 


F  ■+■  Constant  <  1  +  cosift 


(y  +  1)  +  q{y  -  1) 


Applying  L’ Hospital's  rule,  the  above  relation  yields  F-  constant.  There¬ 
fore,  if  k  >  q  as  X  +  0,  the  pressure  approaches  a  constant. 

If  k  -  q,  x  -*■  0  corresponds,  again,  to  <f>  -*■  tt/2;  however  ,  since  n  is 
constant,  F  as  x  -►  0. 

If  k  <  q,  then  k  -  q  <  0  and  n  -*■  “  as  x  -*■  0  since  n  ^  x^  ~  q.  But 
X  *  0  corresponds  to  ♦  -  sin  "hw/fk  *  q(Y  -  1));  therefore,  vl  .. 

From  Equation  34  one  can  see  that  the  distribution  of  F  with  x  as 
x  -v  0  is  determined  to  a  great  extent  by  the  behavior  of  nn‘  as  x  +  0. 

Since  n  ^  x  "  q,  then  on’  *  x2(k  “  q)  “  1  and,  therefore, 

if  k  -  q  >  1/2,  then  nn1  -*  0  as  x  -*•  0; 

if  k  -  q  -  1/2,  then  nn1  =  constant; 

if  k  -  q  <  1/2,  then  nn’  -*•  m  as  x  -*■  0. 

Consider,  first,  k-q  >  1/2  and  examine  the  limit  of  F'  as  x  ^  0.  Noting 

that  z,  ■+  F  one  has 

v'  -  XT  [(?  -  Is)  y’/y] 

applying  L’Hospital's  rule 

P  -  Hr  CC  -  F’)  but  f  =  W'/i  *  CG2nn\)/c 
ar^,  therefore,  as  x  -*■  0, 

-►  P  +  (G^nn1)/^  and  P'  ■+■  y/(l  -y  )  Grin’  -*■  0 

Next  consider  the  limit  of  F'  as  x  0  for  k  -  q  =  i/2.  One  still  has 
S  ■+  P  as  x  -*■  0  and,  therefore, 

p  -►  — 1  —  j^(5»  -  P  )  +  — ■ ~  Constant 

Finally  consider  the  limit  of  F’  for  k  -  q  <  1/2.  As  before  5  -*•  F  as 
x  >  0;  however,  nn’  ->  »,  therefore, 

p’  -  Jx  G2nn'  . 

1  -  Y  P 
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E 


as  x  -*■  0. 

Equation  34  also  indicates  that  no  matter  what  the  relation  of  k  to 
q  that  P  ®  as  c  -*•  ^  or  in  other  words,  the  slope  of  the  curve  for 
V  as  a  function  of  x  is  infinite  as  the  flaw  exit  for  all  choked  flows. 
Figure  5  summarizes  all  of  the  above  results. 


Figure  5.  $  versus  x  and  7  versus  x  for  various  k  and  q  relations 


EFFECT  OF  FRICTION  ON  PRESSURES  PREDICTION  IN  CRACKS  AND  DEBONDS 


Analysis  of  pressure  distributions  in  burning  flaws  have  been  investi¬ 
gated  previously  at  ERDE  ^  ass  lining  that  the  combustion  gases  were  in¬ 
compressible.  They  determined  the  pressure  distribution  assuming  the 
standard  pressure  drop  in  a  duct  due  to  friction.  Although  the  gases 
are  both  viscous  and  compressible,  a  degree  of  satisfaction  was  obtain¬ 
ed  when  the  analysis  was  compared  with  limited  experimental  data  from 
tests  in  which  the  cracks  were  not  choked.  This  gives  rise  to  a  require¬ 
ment  to  evaluate  the  effects  of  both  compressibility  and  friction  to  see 
whether  either  can  generally  be  neglected. 

Figures  6  and  7  show  the  predicted  pressure  distributions  from  the 
analysis  of  the  first  portion  of  this  cliapter  in  both  a  crack  and  debond 
of  similar  geometry  for  a  relatively  high  burning  rate  propellant.  (The 
burning  rate  is  higher  than  for  TP  H-1011.)  The  geometry  was  purposely 
chosen  so  that  in  neither  case  would  the  flow  be  choked  to  the  extent 
that  the  exit  plane  pressure  was  higher  than  the  chamber  pressure.  As 
can  be  seen,  for  both  the  crack  and  the  debond,  a  considerable  error 
would  be  introduced  by  neglecting  the  friction;  however,  the  major  reason 
for  the  large  over  pressures  in  each  flaw  was  due  to  compressibility  ef¬ 
fects  . 

For  the  case  of  the  debond  with  its  lower  over  pressures  and  lower 
mass  flow  the  effect  of  friction  did  not  greatly  affect  the  Mach  number. 

For  the  crack  with  its  higher  mass  flow  and  choked  condition,  the  effects 
of  friction  cause  a  greater  rate  of  change  in  the  Mach  number  near  the 
exit  of  the  crack. 

The  results  illustrated  definitely  support  the  inclusion  of  both 
friction  and  compressibility  in  evaluating  the  fluid  mechanics  of  burning 
in  flaws  in  propellant  grains. 
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DISTANCE  ALONG  DEBOND  ~  INCHES 


Figure  7.  Pressure  and  Mach  number  distributions  along  the  length 

of  a  constantly  diverging  debond  with  and  without  friction, 


MACH  number 


EVALUATION  OF  TWO-DIMENSIONAL  MODELS  OF  FLOW  IN  A  FISSURE  WITH  COMBUSTION 


Presented  in  preceding  sections  was  a  one-dimensiona?  model  for  the 
flow  and  pressure  distribution  in  a  burning  fissure.  It  i  'rtant  to 
determine  in  what  limits  this  one-dimensional  representation  yields  a  re¬ 
liable  representation  of  the  actual  flow.  The  following  analyses  are  in¬ 
tended  to  answer  such  questions.  First,  the  limiting  geometry  is  determined 
in  which  compressibility  effects  become  ;mportant.  Since  flow  from  the 
transpiring  burning  surface  is  rotational,  it  is  also  appropriate  as  a 
second  problem  to  determine  the  effects  of  vorticity  on  flow  in  the  burning 
crack. 

Effects  of  Compressibility 

The  following  is  a  two-dimensional  perturbation  solution  of  the  com¬ 
pressible  flow  field  in  a  sharp  crack  with  arbitrary  mass  addition  at  the 
wall.  A  basic  assunption  is  that  the  Mach  number  of  the  flow  at  the  trans¬ 
piring  wall  is  small;  the  wall  Mach  number  thus  serves  as  a  convenient 
perturbation  parameter.  Another  assumption  is  that  the  flow  is  irrotation- 
al  and  this  severely  affects  the  type  of  boundary  conditions  which  can  be 
satisfied  in  the  solution.  In  a  later  section  the  effects  of  rotationality 
are  assessed.  In  particular,  the  condition  of  normal  effluxe  at  the  wall 
must  be  relaxed  to  secure  an  irrotational  solution.  It  is  readily  demon¬ 
strated  that  an  irrotational  flow  in  geometrical  situations  of  the  present 
type  is  only  possible  with  normal  mass  injection  for  a  special  distribution 
of  wall  sources.  If,  for  example,  the  burning  is  uniform,  then  the  flow 
cannot  be  irrotational  unless  a  velocity  component  parallel  to  the  trans¬ 
piring  wall  is  allowed. 

Figure  8  illustrates  the  assumed  geometry.  Angle  0Q  is  the  half  angle 
of  the  triangular  sharp  crack.  Assuming  irrotational  flow,  the  problem  may 
conveniently  be  represented  by.  a  velocity'  potential  4>  such  that 

v$  =  u  =  ue  +  veQ  (50) 

r  o 

where  u  is  the  flow  velocity  vector.  The  governing  equation  for  steady 
flow  is,  therefore, 
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v2$  (v$*v<j>) 


(51) 


=  v*.v  •  1  j 

where  all  quantities  are  nondimens ionali zed  by  a  characteristic  length 
representative  of  the  length  of  the  crack  and  the  stagnation  speed  of 
sound  of  the  combuaion  gases  aQ.  The  boundary  condition  is  that  the 
normal  components  of  flow  at  the  wall  is 

v  =  i  vf (r)  on  e  =  ±  eQ  (52 

where  v  is  a  Mach  number  representative  of  the  flow  at  the  wall  and  f(r) 
is  the  distribution  of  sources  along  the  wall. 

In  the  plane  polar  coordinates,  the  full  problem  for  compressible 
irrotational  flow  is 
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V  +  ?  *r  *  7  *8e 


. 2  *  .  Vee  . 

*r  V  +  ^7“ 

MS  O0 

_JLi  +  2  *  + 

r3  rz  r 


$0  =  +  v[rf(r)]  on  e  =  +e( 


(S3) 


Hie  foim  of  the  boundary  condition  suggests  a  perturbation  approach  with 
v  as  the  small  parameter.  Put 

+  -  v*(1)  +  vV2)  +  ...  (54) 

and  collect  various  orders  of  v: 

0(v): 

♦S'  +  i  +  4  *QCP  -  0 


-rr  T  r  r  '  7  06 


[  =  +  rf  (r)  on  9  =  + 


(55) 
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r  t 
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and  so  on.  Note  that  the  correction  for  compressibility  enters  at  the 
third-order  term.  The  second-order  solution  is  trivial.  Thus,  as  long 
as  v  is  sufficiently  small,  the  solution  will  correspond  closely  to  the 
incompressible  case. 

Solution  of  the  incompressible  problem  (Equation  55)  is  simple. 

The  eigenfunctions  corresponding  to  the  Laplace  equation  in  plane  polar 
coordinates  are 


=  (A^  sin  k0  +  cos  k0)  (C^r*  +  D^r  K) 

where  <  is  an  integer. 

Solution  to  the  present  problem  may  be  constructed  as  an  infinite 
series  of  these  eigenfunctions: 


CD 

=  53  (A  sin  t<0  +  B  cos  k0)  (C  rK  +  D  r  IC) 


The  coefficients  may  be  determined  by  matching  the  boundary  conditions 
(Equation  55).  Note  by  synmetry,  0,  thus 


0  =  +  0 


=  +  rf(r)  =  +  £  k  sin  k0-  (C  r  +  D  r  K) 
k=1  ,  0  K  * 

o 


(60) 


coefficients  C  ^  D  ^  are  readily  determined  if  f(r)  can  be  expanded  in  a 
Taylor  series.  The  simplest  case  is  for  uniform  transpiration,  f(r)  =  1. 
Then  all  terms  in  the  series  must  vanish  except  for  k  -  1  and  D^  =  0. 

Thus 

r  =  (sin  0  T  C,  r 
v  o  i 

and 

C  =  ^ 

L1  sin  G 

o 

The  velocity  potential  for  this  case  is 
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and  f(r)  =  1  which  shews  that  the  magnitude  of  the  velocity  is  constant. 
+  0(v3) 


u 


sin  6. 


(62) 


and  depends  only  on  the  wall  Mach  number  and  the  crack  half-angle  0Q. 
There  is  a  restriction  on  0Q.  Note  if  0Q  0(v),  then  the  perturbation 
scheme  breaks  down  and  compressibility  effects  are  likely  to  be  important 
For  typical  propellants, 

10"3  <  v  <  10" 2 


so  that  the  present  solution  is  valid  for  cracks  w: th  9q  >  5°.  For 
narrower  cracks,  compressibility  effects  must  be  represented  obviously 
since  the  Mach  number  within  the  crack  approaches  unity.  The  solution 
implies  that  the  flow  may  be  choked  for  a  sufficiently  narrow  crack, 
as  already  indicated  in  the  one -dimensional  solution. 


Effects  of  Vorticity 

It  is  appropriate  to  question  the  validity  of  the  one-dimensional 
analysis  in  representing  the  rotational  flow  generated  by  burning  of  the 
walls  of  a  long  crack.  As  shown  in  the  preceding  analysis  attempts  to 
construct  a  two-dimensional  solution  with  the  assumption  of  irrotational 
flow  do  not  properly  satisfy  the  boundary  conditions.  Such  solutions 
exhibit  a  component  of  velocity  parallel  to  the  surface  of  the  burning 
propellant.  The  burning  process  requires  that  the  velocity  be  normal  to 
the  surface.  Thus  vorticity  is  generated  in  the  burning  process  and  is 
transported  along  the  streamlines;  the  flow  is  rotational. 

Assuming  a  sufficiently  wide  crack,  the  flow  nay  be  assumed  incom¬ 
pressible  as  shown  in  the  preceding  analysis.  Since  the  flow  is  rotation 
al,  it  is  convenient  to  work  with  the  stream  function  ^  such  that 


3  ip 
3X 


u  =  r  30 


3i p 


V  = 


ii  e 

3r  0 


for  cartesian  and  polar  coordinates  respectively. 
The  governing  equation  is 

Vzi|)  = 


(63) 


(64) 


lere  u  is  the  vorticity.  The  vorticity  vector  is  of  course  perpendicular 
to  the  flow  since  planar  flow  is  assumed. 
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Solutions  are  readily  determined  for  the  rectangular  geometry  shown 
in  Figure  9.  The  boundary  value  problem  becomes 


h*k=- 

v  V 

U>(0) 

(65) 

3l p  .  rr  ' 

if  =  1  ^(y. 

on  x  =  +_  a 

(66) 

U  =  o 

3x 

on  x  =  +  a 

(67) 

where  the  flow  is  forced  to  be  normal  to  the  burning  surface  as  required 
by  the  combustion  process.  Solutions  may  be  determined  for  arbitrary 
burning  rate  as  reflected  by  the  function  f(y).  v  is  representative  of 
the  Mach  number  of  combustion  products  near  the  origin  at  the  crack  tip. 
Solutions  of  the  type  required  here  may  be  found  by  assuming  cd  to  be 
proportional  to  Put 

id  =  c2i l>  (68) 

where  c  =  constant. 

This  assumption  is  valid  if  dissipation  within  the  flow  field  is 
neglected.  For  uniform  burning,  f(y)  -  1,  the  appropriate  solution  is 
simply 


*  =  -vy  sin(^)  (69) 

which  assumes  no  burning  at  the  crack  tip.  Streamlines  for  this  solution 
are  shown  in  Figure  9. 
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III.  PARAMETRIC  STUDY  OF  THE  PRESSURE  DISTRIBUTIONS 
FOR  BURNING  CRACKS  AND  DEBONDS 

Presented  herein  is  a  parametric  study  of  the  influence  of  geomet¬ 
ric  and  propellant  characteristics  on  the  pressure  distribution  in  cracks 
and  debonds.  The  study  is  carried  out  primarily  for  TPH-1011  propellant; 
however,  the  effects  of  higher  burning  rates  are  indicated. 

EFFECT  OF  FLAW  LENGTH  AND  EXIT  AREA  ON  CRACK  TIP  PRESSURE 

Although  the  entire  distribution  of  pressure  along  the  length  of  a 
flaw  is  necessary  to  calculate  the  surface  energy  of  the  flaw,  the  pres¬ 
sure  most  characteristic  of  the  loading  is  the  tip  pressure.  This  is 
the  pressure  at  the  location  of  propagation  and  the  location  of  the  maxi¬ 
mum  burning  rate  as  well.  Therefore,  to  illustrate  the  effects  of  flaw 
geometry  on  the  pressures  within  a  burning  flaw,  the  flaw  tip  pressure 
is  presented  as  a  function  of  flaw  length.  Assuming  a  triangular  or 
uniformly  diverging  flaw  for  illustrative  purposes,  the  exit  crack  width 
serves  to  illustrate  the  effects  of  exit  area. 

Figures  10  and  11  indicate  the  effects  of  flaw  length  and  exit  width 
for  a  crack  and  a  debond  respectively.  The  propellant  properties  used 
are  those  of  TPH-1011.  It  is  assumed  that  the  main  chamber  pressure  is 
1000  psia  and  that  the  crack  tip  width  is  0.01  inches. 

For  both  cases  it  is  seen  that  the  flaw  tip  pressure  increases 
rapidly  with  increasing  length  and,  further,  that  the  pressure  increases 
with  decreasing  angle  of  divergence  (decreasing  flaw  width) . 

EFFECT  Or  ASSUMED  FLAW  TIP  DIMENSION^ 

Numerical  solutions  of  the  governing  one-dimensional  equations  pre¬ 
sented  in  Section  I I -A  for  cracks  and  debonds  cannot  be  obtained  if  the 
crack  is  assumed  sharp.  In  addition,  there  is  no  such  thing  as  a  truly 
sharp  crack.  However,  it  would  appear  reasonable  to  evaluate  the  ef¬ 
fect  of  the  assumed  tip  geometry  numerically  as  it  approaches  zero  to 
ascertain  the  form  of  the  singularity  at  zero  tip  area. 
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Crack  Length  'v  Inches 

Figure  10  Crack  tip  pressure  versus  length  for  ieveral  Crack  geometries 
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Figures  12  and  13  present  typical  results  for  both  a  triangular 
crack  and  a  triangular  debond  of  otherwise  constant  geometry.  The 
geometry  chosen  was  a  length  of  two  inches  and  an  exit  width  of  0.05  inches . 
The  propellant  was  assumed  to  be  TPH-1011  and  the  chamber  pressure  1,000  psia. 
For  both  types  of  flaws  there  were  small  changes  in  tip  pressure  for  a 
decrease  in  tip  width  from  0.01  inches  to  0.001  inches.  Further  decreases 
in  tip  width  down  to  10"^  inches  caused  changes  of  only  400  psi  for  the 
case  of  the  crack  (Figure  13}  and  200  psi  for  the  debond  (Figure  12) . 

EFFECT  OF  BURNING  RATE  LAW 

Figures  14  and  15  indicate  the  effects  of  propellant  burning  rate 
coefficients,  GDC,  on  the  flaw  tip  pressure  over  a  range  of  flaw  lengths. 

The  analysis  assumed  constant  tip  and  exit  plane  widths  and  a  constant 
chamber  pressure.  C  =  0.0047  is  equivalent  to  TPH-1011.  As  can  be  seen, 
the  propellant  burning  rate  has  a  strong  effect  on  the  resulting  pressures. 
While  the  lower  burning  rate  propellant,  C  =  0.00027  shows  negligible 
increase  in  pressure,  TPH-1011  and  faster  burning  propellants  rapidly 
reach  high  over  pressures.  Thus,  all  conclusions  reached  here  for  TPH-1011 
flaw  tip  pressure  would  be  conservative  estimates  of  what  would  occur  in 
faster  burning  propellants  such  as  the  double  base  propellants  utilized  in 
some  operational  systems. 

Figures  16  and  17  present  plots  of  exit  plane  Mach  lumbers  for  the 
propellants  of  Figures  14  and  15.  The  higher  burning  rate  propellant 
reaches  choked  conditions  for  flaw  lengths  considerably  less  than  those 
considered. 

TYPICAL  INCREASE  IN  FLAW  LBJGTH  WITH  TIME,  ASSUMING  NO  MECHANICAL  PROPAGATION 

The  computer  program  described  in  Appendix  B  will  account  for  all 
changes  in  flaw  geometry  as  a  function  of  time,  assuming  that  the  flaw 
does  not  propagate  mechanically.  Typical  results  for  cracks  and  debonds 
are  presented  in  Figures  18  and  19  respectively.  As  can  be  seen,  for  the 
constant  initial  width  flaws  assuming  TPH-1011  propellant  with  a  chamber 
pressure  of  1,000  psia,  the  rate  of  increase  in  length  with  time  varies 
directly  with  initial  crack  length.  This  is  as  would  be  expected  since 
it  correlates  directly  with  the  pressure  loadings.  The  crack  growth 
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CRACK  LENGTH  ~  INCHES 


Figure  16.  Exit  Mach  number  versus  crack  length  for  several  propellants 
at  constant  tip  and  exit  width 


Exit  Mach  Number 


Figure  17.  Exit  Mach  number  versus  debond  length  for  several  propellants 
at  constant  tip  and  exit  width. 
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Figure  18.  Increase  in  crack  length  versus  time  for  several  initial  lengths. 
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TIME  ~  SECONDS 


due  to  burning  alone  can  safely  be  said  to  be  sufficiently  slow  so  that 
it  offers  no  great  problem. 

HALF-TIME  FOR  TIP  PRESSURE  VERSUS  BURNING  RATE  COEFFICIENT 

A  characteristic  of  all  burning  flaws  is  that  the  pressure  decreases 
rapidly  with  time  due  to  the  geometrical  changes  occurring.  To  illustrate 
this  effect,  Figures  20  and  21  indicate  the  time  required  to  reduce  the 
flaw  tip  pressure  by  a  factor  of  two  as  a  function  of  the  coefficient, 

C,  in  the  burning  rate  law.  The  burning  rate  law  is  that  for  TPH-1011  with 
coefficient  changed.  As  would  be  expected,  the  half  life  drops  off  at 
a  faster  rate  for  the  crack  than  the  debond  due  to  the  fact  that  it  is 
burning  on  both  surfaces.  Although  extremely  high  pressures  may  exist 
initially  in  a  flaw,  the  defect  geometry  changes  at  a  near  constant  rate 
due  to  the  small  exponent  in  the  burning  rate  law  as  can  be  seen  in 
Figures  18  and  19. 

EFFECT  OF  VARIABLE  GEOMETRY.  CIRCULAR  DFBQNDS 

Figures  22  and  23  represent  the  types  of  situations  which  might 
result  in  a  head  end  debond  around  a  rocket  motor  igniter.  As  can  be 
seen,  the  smaller  the  igniter  diameter,  the  greater  the  difficulty  for 
equally  deep  debcnds. 


Time  for  Tip  Pressure  to  Reduce  to  One  Half  Initia 
Pressure  ~  sec 
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Figure  20.  Time  for  tip  pressure  to  reduce  to  half  initial 

pressure  versus  burning  rate  coefficient  C  for  a  given  initial 
crack  geometry. 
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Figure  21.  Time  for  tip  pressure  to  reduce  to  half  initial 

pressure  versus  burning  rate  coefficient  C  for  a  given 
initial  debond  geometry 
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ach  Number  P  _  -  P  ,  nsi 


Length  Along  Circular  Debond  ~  inches 

Figu:  ?  23.  Pressure  and  Mach  number  distribution  along 
circular  debond 
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IV.  TIME  HISTORY  OF  FRACTURE  OR  DEBONDING  IN 


VISCOELASTIC  MATERIALS 


Stated  in  perhaps  its  simplest  form, the  problem  to  be  solved  from 
the  solid  mechanics  point  of  view  is  the  following:  given  the  geometric 
configuration  of  the  crack  or  debond  in  the  rocket  grain  and  the  pressure 
distribution  inside  this  crack  or  debond,  both  as  known  functions  of 
time,  determine  the  instant  at  which  fracture  is  initiated  and,  at  least 
for  small  times,  the  time-history  of  the  fracture.  Stated  in  such  a 
maimer,  the  problem  may  appear  to  be  a  relatively  simple  one.  Such, 
unfortunately,  is  not  the  case  due  primarily  to  the  difficulties  en¬ 
countered  in  mathematically  analyzing  a  body  with  a  general  cracked 
geometry.  When  the  geometry  changes  with  time  as  in  this  problem 
once  fracture  has  been  initiated,  or  the  loads  are  applied  to  the  grain 
dynamically,  inertia  effects  may  become  of  sufficient  importance  so  that 
they  can  no  longer  be  neglected.  Thus,  a  new  dimension  of  difficulty  is 
added  to  the  problem.  Finally,  it  should  be  noted  that  the  materials 
one  must  work  with  when  investigating  crack  propagation  in  solid  oro- 
pellant  rocket  motors  must  be  characterized,  at  best,  as  being  linearly 
viscoelastic.  Unlike  elastic  materials  in  which  the  stress  is  a 
function  of  strain  alone,  the  rate  dependency  of  viscoelastic  materials 
means  that  the  stress  depends  upon  the  entire  time-history  of  the  strain. 
This,  of  course,  complicates  the  analysis  still  further. 

Fortunately,  however,  the  continuum  theory  of  fracture  has  received 
a  great  deal  of  attention  in  this  century  since  its  beginnings  in  the 
classical  work  of  Griffith  Cl 5)  and  much  progress  has  been  made  due 
in  large  part  to  the  improvement  of  mathematical  techniques  and  the 
advent  of  the  high  speed  computer.  This  progress  has  been  largely  con¬ 
centrated  in  the  areas  of  fracture  of  brittle  elastic  materials  and 
more  recentl> ,  of  fracture  of  ductile  elastic  materials  where  the  effects 
of  plastic  deformation  in  a  region  surrounding  the  crack  tip  are  con¬ 
sidered.  However  with  regard  to  the  present  study,  it  should  be 
mentioned  that  a  great  deal  of  work  has  been  done  recently  on  developing 
methods  of  analysis  for  viscoelastic  structures  of  uncracked  geometries. 
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In  addition,  many  of  the  concepts  of  fracture  mechanics  develoi*eu  for 
use  with  elastic  materials  are  still  valid  for  use  with  viscoelastic 
materials. 

Because  of  the  noted  lack  of  information  available  specifically  on  visco¬ 
elastic  fracture,  an  approach  to  the  solid  mechanics  part  of  this  program  has 
been  chosen  which  is  based  upon  the  development  and  analysis  of  relatively 
simple  models  capable  of  describing  the  important  aspects  of  viscoelastic 
fracture  and  debonding  at  least  insofar  as  the  viscoelastic  dependence 
is  concerned.  The  information  obtained  in  analyzing  these  models  could 
then  be  used  to  provide  qualitative  results  in  the  prediction  of  fracture. 
With  these  comments  in  mind,  let  us  consider  the  progress  that  has  been 
made,  first,  in  the  field  of  cohesive  fracture  and  then  in  the  closely 
related  (from  a  continuum  mechanics  point  of  view)  field  of  adhesive 
debonding . 

The  first  work  to  be  mentioned  is  the  investigation  by  Swanson 
of  fracture  in  a  linearly  viscoelastic  tubular  rocket  propellant  grain. 

In  this  study,  the  cylindrical  grain  was  assumed  to  be  of  infinite 
length  and  to  have  a  radial  crack  at  its  centerbore  running  the  full 
length  of  the  cylinder.  Time -dependent  critical  stress  intensity 
factors  K^(t)  were  then  used  as  criticality  conditions  for  the  cylinder 
subjected  to  a  time -dependent  internal  pressure  loading. 

The  critical  stress  intensity  factors  K^(t)  to  be  used  were 
obtained  from  laboratory  tests  on  specimens  of  Hercules  Incorporated 
designated  EJC  solid  propellant.  These  specimens  were  in  the  form  of 
solid  right  circular  cylinders  one  inch  in  diameter  by  three  inches 
long  with  a  crack  three -sixteenths  of  an  inch  deep  machined  circum¬ 
ferentially  around  the  specimen  by  knife  blade  on  a  lathe.  Bonding 
these  specimens  to  rigid  end  plates,  they  were  then  subjected  to 
several  constant  rate  tensile  tests  at  various  cross-head  speeds.  In 

addition,  a  series  of  tests  was  also  conducted  in  which  the  crosshead  speed 
was  changed  from  one  constant  value  to  another  during  the  course  of  the 
test.  In  each  of  these  tests  the  time  to  failure  and  the  failure  load 
were  recorded.  From  the  load  at  failure,  it  was  a  simple  matter  to 

calculate  the  stress  intensity  factor  at  failure  from  the  analysis 

r  j.7i 

given  by  Bueckner  ^  J  for  an  elastic  cracked  cylindrical  specimen. 

Bueckner's  results  can  be  summarized  as 
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where  6ngt  is  the  axial  force  divided  by  the  net  area  of  the  notched 
section,  D  is  the  outside  diameter  of  the  specimen,  d  is  the  diameter 
of  the  notched  section  and  F(d/D)  is  a  numerical  factor  listed  in 
Reference  17.  It  is  clear  from  Equation  70  and  the  elastic- visco¬ 
elastic  correspondence  principle  that  this  stress  intensity  factor  applies 
to  the  viscoelastic  test  specimens  as  well.  Hence,  plotting  at 
failure,  i.e.  K^,  versus  the  time  to  failure  for  the  various  tests 
performed,  a  plot  of  K^(t)  was  obtained  for  the  EJC  propellant.  This 
is  shown  in  Figure  24. 

Now  if  a  viscoelastic  analysis  of  a  particular  cracked  structure  of 
EJC  propellant  can  be  performed  so  that  the  time -dependent  stress 
intensity  factor  K^t)  can  be  found,  the  time  to  fracture  is  given  by 
the  intersection  of  the  Kic(t)  and  K^t)  curves.  The  difficulty  here 
lies  in  the  fact  that  since  K^(t)  was  plotted  for  constant  loading 
rates,  our  predictions  of  time  to  fracture  are  necessarily  restricted 
to  structures  undergoing  at  least  an  approximately  constant  loading 
rate.  Fortunately,  however,  use  of  the  spherical  flaw  model  (to  be 
discussed  later  in  this  section)  and  laboratory  tests  have  shown  the 
time-dependent  critical  stress  intensity  factor  approach  to  be  relatively 
insensitive  to  a  wide  range  of  fast-slow  variable  loading  conditions. 

Using  a  finite-element  stress  analysis  method,  a  plane-strain 
elastic  analysis  was  made  for  a  radial  crack  at  the  centerbore  of  a 
tabular  rocket  propellant  grain.  The  analysis  was  made  to  calculate 
the  strain  energy  in  the  propellant  for  a  crack  depth  one-half  panel 
greater  and  one-half  panel  less  than  the  actual  crack  depth  under 
consideration.  Having  calculated  each  of  these  strain  energies,  the 
strain  energy  release  rate  with  respect  to  crack  surface  area  is 
then  calculated  from  a  finite  difference  approximation  to  the  formula 
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Plot  of  time -dependent  critical  stress 
intensity  factor  K,  (t). 


where  U  is  the  strain  energy  and  A  is  the  crack  surface  area.  Having 
calculated  this  expression,  one  can  finally  obtain  the  stress  intensity 


factor  Kj  as 


EG, 


1-v 


(72) 


where  E  is  Young's  modulus  and  v  is  Poisson's  ratio. 

Once  having  obtained  the  stress  intensity  factor  for  the  elastic 
problem,  it  is  possible  to  solve  for  the  time -dependent  stress  intensity 
factor  K1(t)  for  the  corresponding  viscoelastic  problem  through  the 
use  of  Schapery's  quasi-elastic  approximation  In  fact,  this  results 

in 

K‘2(t}  ‘  (K>2)  (¥9  C73) 

where  Erel(t)  is  the  stress  relaxation  modulus,  if  the  internal  pressure 
is  applied  as  a  step  function  in  time.  When  the  pressure  does  not 
vary  with  time  in  this  manner  but  in  some  more  general  way,  the  fact 
that  we  are  considering  a  linear  viscoelastic  material  and  that  Kj  is  a 
linear  function  of  the  applied  pressure  loading,  allows  one  to  use 
Equation  73  with  the  Duhamel  integral  to  obtain  KjCt)  for  this  more 
general  time-varying  pressure  loading.  Plotting  the  resulting  Kj(t) 
for  a  particular  time-varying  internal  pressure  on  the  same  graph  with 
the  K^(t)  curve  obtained  from  the  laboratory  tests  on  cracked  specimens, 
the  time  to  failure  is  predicted  as  the  point  of  intersection  of  the 
two  curves. 

To  verify  the  analytical  work,  a  high-rate  hydrotest  was  conducted 
on  a  structural  test  vehicle  (STV)  containing  a  radial  crack.  Crack 
propagation  was  produced  during  the  test  and  the  time  of  its  initiation 
was  obtained  from  instrumentation  recordings  of  grain  deflections  and  the 
pressure  rate.  The  measured  pressure-time  curve  was  then  used  in  the 
previously  mentioned  finite -element  analysis  of  a  viscoelastic  grain 
to  obtain  the  time-dependent  stress- intensity  factor  for  the  STV,  This 
calculated  K^t)  curve  was  plotted  in  Figure  25  along  with  the 
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experimentally  obtained  Kic(t)  curve.  The  irregularity  in  the  K^t) 
curve  was  believed  to  be  caused  by  entrapped  air  in  the  hydraulic 
system.  The  intersection  of  the  two  curves  in  Figure  25  is  the  predicted 
time  to  failure  and  is  seen  to  he  within  five  per  cent  of  the  observed 
time  to  failure. 

Although  in  the  one  test  conducted  the  stress  intensity  factor 
approach  outlined  above  was  quite  successful  in  predicting  the  time  to 
failure,  two  limitations  of  this  method  do  arise.  The  first  is  that  a 
uniform  pressure  distribution  was  assumed  to  act  within  the  crack  whereas 
the  analysis  of  Section  III  shows  that  the  actual  pressure  distribution 
in  a  burning  crack  is  non-uniform.  The  effects  of  such  a  non-uniform 
pressure  distribution  on  crack  instability  can  be  estimated,  however, 
by  the  results  of  Appendix  C,  thus  causing  no  substantial  difficulties. 

A  more  serious  objection  is  that  the  critical  stress  intensity  factor 
approach  is  not  applicable  to  the  second  part  of  our  problem;  that  is, 
the  prediction  of  initial  velocities  of  propagation.  Hence,  we  have 
considered  a  more. general  approach  to  the  problem  based  on  the  thermo¬ 
dynamic  power  balance.  This  balance  can  be  written  as 

I  =  F  +  -  2D  +  SE  +  K  (74) 


where  I  is  the  power  input  of  the  applied  loading  at  the  boundaries  of 
the  system,  r  is  the  rate  of  increase  of  the  free  (strain)  energy, 

2D  is  the  dissipation,  SE  is  the  rate  of  increase  of  the  surface  energy 
and  K  is  the  rate  of  increase  of  the  kinetic  energy.  In  principle, 
this  power  balance  can  be  applied  directly  to  any  crack  configuration, 
loading  and  material  of  interest.  In  practice,  however,  computational 
difficulties  are  encountered  in  the  application  of  the  power  balance 
to  realistic  geometries  and  materials.  This  has  led  us  to  consider 
William’s  model  •  of  a  spherical  flaw  which  incorporates  a  simplified 
crack  geometry  and  the  power  balance  Equation  74  tc  predict  both  the 
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Figure  25.  Comparison  of  predicted  and  observed  crack 
propagation  in  pressure  test  of  STV. 


i 


time  to  fracture  and  its  initial  velocity  of  propagation.  Aside  from 
the  relative  ease  with  which  one  can  calculate  the  various  terms  of  the 
power  balance  (Equation  74)  for  this  geometry,  further  motivation  for 
studying  the  spherical  flaw  was  provided  by  the  fact  that  the  criticality 
condition  for  the  elastic  case  is  extremely  similar  to  that  found  for 
the  more  realistic  crack  geometries  of  Griffith^15'*  or  Sneddon^20^. 

Since  this  model  has  been  extensively  investigated  in  the  literature, 
it  is  sufficient  here  only  to  outline  its  important  features  as  regards 
this  study. 

In  his  original  work  on  the  spherical  flaw,  Williams^'^  neglected 
inertia  effects  and  considered  the  surface  of  the  flaw  to  he  stress 
free  while  subjecting  the  outer  boundary  of  the  hollow  sphere  to  four 
typical  inputs,  namely,  constant  stress  or  displacement  and  constant 
stress  or  displacement  rate.  For  each  of  these  four  loadings  he  was 
able  to  arrive  at  an  expiession  from  which,  to  calculate  the  time  to 
fracture.  However,  for  this  simple  model  difficulties  were  encountered 

4 

in  solving  the  nonlinear  integral  equation  for  the  time  history  of  the 
flaw  growth.  For  example,  an  initial  velocity  of  propagation  was  found 
for  the  car .  where  the  outer  sphere  boundary  was  subjected  to  a  constant 
displacement  rate,  but  in  the  limiting  case  of  instantaneous  propagation, 
that  is,  time  to  fracture  t^  =  0,  this  initial  velocity  was  discovered 
to  be  infinite.  Tins  rather  unrealistic  result  was  attributed  to  the 
fact  that  inertia  effects  were  neglected  in  the  analysis.  Although  this 
omission  may  cause  little  error  in  a  highly  viscous  material  subjected 
to  slow  loading  rates,  this  will  not  be  the  case  in  more  elastic  materials 
or  when  the  applied  loading  rates  become  appreciable.  This  was,  in  fact, 
borne  out  by  a  subsequent  analyses of  the  spherical  flaw  in  which  the 
kinetic  energy  contribution  to  the  power  balance  was  included.  Here  it 
was  shown  that,  for  the  case  of  a  constant  displacement  rate  loading 
at  the  sphere's  outer  boundary,  the  power  balance  fracture  criterion 
predicted  a  zero  initial  velocity  of  propagation  of  the  flaw  when 
fracture  was  initiated  instantaneously.  Thus,  the  effect  of  including 
the  kinetic  energy  in  this  particular  example  was  to  exhibit  a  smooth 
transition  during  the  acceleration  from  zero  initial  flaw  velocity. 
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Because  of  its  wide  applicability  and  its  ability  to  predict  the  time- 
history  of  fracture,  the  energy  balance  will  be  applied  as  the  criterion 
for  fracture  in  our  further  studies.  Thus,  the  spherical  flaw  has 
served  several  useful  purposes  toward  this  end.  Namely,  it  has  provided 
us  with  a  better  understanding  of  the  subtleties  of  applying  the  power 
balance;  it  has  given  us  some  qualitative  insight  into  the  behavior  of 
crack  propagation  in  linearly  viscoelastic  materials;  it  has  also  provided 
(see  Section  V)  several  useful  estimates  regarding  times  to  failure  and 
the  magnitudes  of  initial  crack  propagation  velocities  in  burning  rocket 
grains,  All  of  these  will  be  of  benefit  as  more  sophisticated 
models  are  developed  and  solved  in  the  process  of  obtaining  the  final 
desired  quantitative  results. 

Unlike  the  field  of  cohesive  fracture,  little  work  has  been  done 

in  the  field  of  adhesive  fracture  (debonding)  from  a  continuum  mechanics 

approach,  even  for  the  case  of  bonded  elastic  materials.  In  fact,  it 

was  only  recently  pointed  out  that  since  real  adhesive  interfaces,  like 

* 

real  materials,  contain  small  cracks  (debonds)  which  give  rise  to 
stress  concentrations,  the  Griffith  approach  to  cohesive  fracture 
could  then  be  extended  to  the  study  of  adhesive  debonding.  Williams . 
for  example,  showed  that  when  using  the  power  balance  as  a  criterion 
for  crack  propagation,  the  only  difference  mathematically  between  the 
phenomena  of  cohesive  fracture  and  adhesive  debonding  is  in  the 
interpretation  of  the  energy  required  to  create  new  free  (cohesive  or 
adhesive)  surface  area. 

To  this  date,  continuum  mechanics  studies  of  adhesive  failure  have 
centered  almost  entirely  on  studies  of  the  mathematical  singularities 
which  occur  at  the  tip  of  cracks  along  an  interface  of  two  dissimilar 
media  (23-24)  Qr  Qn  methods  for  the  experimental  determination  of  the 
energy  required  to  create  new  free  adhesive  surface  area,  Y  (25-26) ^ 

a 

These  two  points  are  worthy  of  further  discussion. 

As  for  the  case  of  cracks  in  a  homogeneous  solid,  the  stress 
singularities  at  the  tip  of  a  crack  tilong  the  interface  of  two 
dissimilar  media  can  be  found  by  the  linear  theory  of  elasticity. 

However,  in  the  general  case  of  arbitrary  material  constants  the 
singularities  are  not  solely  of  the  square  root  type  as  they  are  for 
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cracks  in  homogeneous  media,  hut  arc  found  to  oscillate  rapidly  near  the 
crack  tip.  For  example,  the  stresses  are  of  the  form^27^ 


■1  [aq  cos (6  tn-S-J  -  B0  sin(3 
±  [A0  sin(0  in^)  +  Bq  cos  (3  in^)] 


+  0(1) 


+  0(1) 


(75a) 


(75b) 


for  the  case  of  a  oond  line  subjected  to  both  tension  and  in-plane  shear. 
In  Equation  75,  3  is  defined  by 


3 


>*!  +  v2&  '  4vj) 
VZ  +  Ml<*3  ”4  v2^ 


(76) 


Aq  and  Bq  are  stress  intensity  factors,  r  is  the  distance  from  the 

crack  tip  along  the  interface,  a  is  the  crack  half  length,  m  and  v  are 

the  shear  modulus  and  Poisson's  ratio,  respectively,  and  the  subscripts 

one  and  two  refer  to  the  different  materials  on  either  side  of  the  debond. 

Note  that  for  the  special  case  where  the  debond  occurs  at  the  interface 

between  a  rigid  material  and  an  incompressible  one,  that  is  u2  r  m  and 

=  0.5,  no  oscillations  occur.  Because  of  the  trig-log  behavior  of  the 

stresses  in  the  general  material  case,  one  cannot  rigorously  compute 

the  strain  energy  directly.  However,  as  a  practical  matter,  the 

oscillating  singularity  is  disregarded 

The  experimental  determination  of  the  energy  required  to  create 

new  free  adhesive  surface  area,  y  ,  has  been  the  subject  of  much 

discussion  and  several  methods  have  been  proposed.  Perhaps  the  best 

of  these  methods,  however,  from  the  combined  viewpoint  of  experimental 

r28) 

convenience  and  ease  of  mathematical  modeling,  is  Williams'1'  '  modifica¬ 
tion  of  a  method  originally  proposed  by  Dannenberg  In  this 

"pressurized  blister  test  "  a  thin  disk  layer  of  soft  material  is  cast 
and  cured  upon  a  relatively  rigid  base  plate  except  for  a  central 
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circular  portion  which  is  prevented  from  bonding.  Pressure  is  then 
introduced  into  this  unbonded  region  causing  the  thin  disk  layer  to 
lift  up  off  the  base  plate  in  a  fashion  similar  to  a  blister.  The 
radius  of  the  unbonded  region  remains  constant  as  the  pressure  is 
slowly  increased  until  a  certain  critical  pressure  level  is  reached. 

At  this  time  the  radius  enlarges,  signifying  an  adhesive  failure 
(debonding)  along  the  interface  of  the  disk  and  base  plate.  Modeling 
the  system  mathematically  and  applying  the  energy  balance  results 
in  a  relftionship  between  the  material  properties  of  the  disk,  the 
critical  pressure  and  y„.  This  result,  combined  with  the  experimental 

a 

data,  ’lews  one  to  arrive  at  a  value  for  y  . 

Since  the  experimentally  determined  values  of  y„  are  critical  to 

a 

the  successful  use  of  the  power  balance  in  predicting  adhesive  debonding, 
it  is  imperative  that  these  values  be  as  accurate  as  possible.  This 
accuracy,  however,  aside  from  laboratory  techniques,  is  heavily  dependent 
15) on  how  well  the  mathematical  model  of  the  blister  test  can  be  made 
to  conform  to  the  actual  experimental  setup.  Thus,  some  effort 
lias  been  expended  in  developing  more  sophisticated  mathematical 
models  of  the  blister  test  to  replace  the  original  analysis.  In  this 
first  analysis  ^22^  the  thin,  flexible  disk  layer  was  modeled  by  linear 
elastic  plate  theory  while  the  base  plate  was  considered  to  be  rigid. 
Since  then  two  more  sophisticated  models  have  been  made.  The  first^2^ 
used  non-linear,  large  deflection  plate  theory  to  evaluate  the  effects 
of  large  deflections  of  the  disk  layer  while  the  second'-'^  was 
concerned  with  the  case  of  two  materials  bonded  together  by  a  third 
rather  than  one  being  cast  directly  upon  the  other.  Let  us  investigate 
each  of  these  models  further. 

It  can  easily  be  shown  using  linear  elastic  plate  theory  and  the 
power  balance  equation  (again  neglecting  inertia  effects)  that  the 
energy  required  to  create  new  free  adhesive  surface  area  in  a  blister 
test  is  given  by^2^ 


^crw0  "  “Ya 


(77) 
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where  pcr  is  the  pressure  needed  to  initiate  the  debonding  and 
is  the  center  deflection  of  the  blister  at  the  time  debonding  starts. 
Thus,  values  of  y  can  be  obtained  from  laboratory  blister  tests  by 
simply  recording  the  pressure  and  central  deflection  of  the  disk  at  the 
time  the  debond  begins  propagating.  A  plot  of  pcr  versus  wQ  is  shown 
in  Figure  26  for  a  test  run  on  polyurethane  cast  and  cured  onto  a 
polished  glass  plate.  It  is  noted  in  this  plot  that  for  small  center 
deflections  the  fracture  initiation  points  fall  quite  closely  along 
a  hyperbola  of  parametric  value  2y  .  However,  for  larger  deflections 
this  is  no  longer  the  case  due  to  the  inability  of  linear  plate  theory 
to  account  for  the  midplane  stretching  in  the  polyurethane  specimen. 
Performing  a  membrane  analysis  and  again  using  the  power  balance,  we  find 
that  if  the  membrane  stresses  are  such  that  they  overshadow  the  bending 
stresses,  yo  is  given  as 

a 


PcrK0 


=  2.4y 


(78) 


Thus  for  larger  center  deflections,  the  plot  of  fracture  initiation 

points  in  Figure  26  should  fall  on  a  hyperbola  of  parametric  value  2.4ya. 

This  is  seen  to  be  the  case.  In  the  transition  between  linear  plate 

(31) 

theory  and  membrane  theory,  an  approximate  solution  by  Berger v  '  can 
be  used  to  include  the  effects  of  both  plate  bending  and  stretching. 

In  this  case,  the  us^  of  the  power  balance  results  in  a  complicated 
expression  from  which  to  find  Ya.  This  expression  must  be  solved 
numerically. 

The  effects  on  the  experimentally  determined  values  of  Y  .due  to 

a 

the  thickness  of  two  materials  being  bonded  together  by  a  third  have 
been  investigated*-"5^  by  approximating  the  interlayer  as  a  Winkler-type 
elastic  foundation.  Using  for  the  analysis  the  standaid  blister  test 
configuration  with  the  addition  of  the  interlayer,  an  application  of 
the  power  balance  results  in  an  expression  for  y  which  can  be  written 
in  the  form 


p  a 

i  cr 

128D 


[1 


(79) 
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Figure  26.  Plot  of  critical  pressure  versus  center  deflection 
from  a  blister  test  of  polyurethane  cast  and  cured 
on  a  glass  plate. 
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Figure  27.  Configuration  of  end  loaded  rocket  propellant  grain 
debonding  axisyninetrically  from  rigid  casing. 


where  pcr  is  the  pressure  at  the  instant  of  fracture  initiation,  a  is 
the  radius  of  the  debond  before  fracture,  D  is  the  flexural  rigidity  of 
the  plate  and 

•[•  v  .IN 

(1  -  2V1) (1  +  v1)]  ^4Dn  J 

In  Equation  80,  v*  and  E*  are  Poisson's  ratio  and  Young's  modulus, 
respectively,  for  the  interlayer  and  h1  is  the  thickness  of  the  inter¬ 
layer.  It  is  clear  that  in  the  limiting  case  of  a  vanishing  interlayer, 
i,e.  h^  -*•  0,  we  have  A  Values  of  the  function  f(A  )  are  given  in 

Table  I,  An  inspection  of  these  values  (in  particular  A  shows 

that  f (A  )  is  the  correction  to  y  due  to  the  infli  nee  of  the  interlayer 

3  3 

thickness. 

It  should  be  mentioned  here  that  both  y  and  y  ,  that  is,  the  energy 
necessary  to  create  new  cohesive  or  adhesive  fracture  surface,  respectively, 
have  usually  been  assumed  to  be  time -independent.  This  assumption  was 
imposed  primarily  as  a  matter  of  analytic  convenience  although  it  has  beer 
widely  recognized  that  this  is  not  the  case  t33-34)^  -j^e  time  dependence 
of  y  and  the  slow  growth  of  adhesive  fracture  observed  in  the  blister 
tests  that  have  so  far  been  made,  suggest  that  time -dependent  dissipative 
mechanisms  take  place  duru.ig  fracture.  Thus,  a  further  refinement  in 
the  mathematical  model  of  *■’. ;  blister  test  is  still  necessary  in  that 
the  thin,  flexible  disk  layer  must  be  modeled  as  being  viscoelastic. 

Work  is  proceeding  in  this  direction. 

Assuming  the  same  materials  and  surface  preparation  as  were  used  in 
the  experimental  determination  of  >  ,  one  can  use  this  now  known  value 
of  y  to  predict  adhesive  debonding  in  a  specimen  of  different  geometry 
under  different  loading  conditions  if  one  can  solve  this  problem  for 
the  potential  energy  release  rate  with  further  debonding.  For  example, 
a  rather  idealized  problem  of  adhesive  debonding  in  a  case-bonded  solid 
propellant  rocket  motor  has  been  investigated  as  an  illustration  of 
the  approach  to  be  used.  The  idealization  considered  is  that  of  a 
finite  length,  elastic  solid  propellant  grain  debonding  in  an  axially 
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symmetric  manner  from  one  end  of  its  enclosing  rigid  case.  The  end 
opposite  the  finite  length  debond  is  considered  to  be  fixed  and  the  grain 
is  assumed  to  be  solid  (100  per  cent  web  fraction  grain).  The  end  of  the 
grain  where  the  debond  takes  place  is  subjected  to  a  normal  pressure  which 
also  acts  in  the  debonds  and  the  dimensions  and  material  properties  are  as 
shown  in  Figure  27.  If  we  further  assume  that  the  kinetic  energy  of  the 
system  cai  be  regarded  as  being  negligible,  the  power  balance  relation 
can  be  written  using  Clapeyron’s  theorem  as 


»U 


(81) 


where  U  is  the  strain  energy  of  the  system  and  A  is  the  surface  area  of 
the  debond.  Thus,  the  problem  becomes  one  of  finding  the  variation  of 
the  strain  energy  with  respect  to  changes  in  the  debond  area.  For  the 
geometry  shown,  this  has  been  accomplished  using  the  finite  element 
method  for  each  of  several  debond  lengths.  The  results  of  the  analysis 
are  shown  in  Figure  28  where  the  strain  energy  is  plotted  against  the 
debond  area.  The  inverse  square  root  of  the  slope  of  this  curve  was 
then  used  with  Equation  81  to  produce  the  parametric  design  curve 
shorn  in  Figure  29.  This  curve  represents  the  debonding  failure 
criterion  for  end  pressure  loadings  on  an  axisymmetric  elastic  rocket 
grain  banded  to  a  rigid  casing. 

Since  it  is  well-known  that  most  solid  rocket  fuels  are  at  best 
linearly  viscoelastic,  the  above  analysis  would  be  improved  considerably 
by  assuming  the  grain  to  be  a  linear  viscoelastic  material.  Although 
it  is  again  theoretically  possible  to  calculate  the  quantities  necessary 
for  use  with  the  power  balance,  the  viscoelastic  nature  of  the  grain 
causes  the  stress  and  strain  fields  to  be  time -dependent,  thereby 
complicating  the  analysis  considerably.  In  this  case  Equation  81 
is  no  longer  applicable.  Rather  we  must  employ  the  more  general  version 
of  the  power  balance  as  given  by  Equation  74.  If  the  various  terms  of 
this  equation  could  then  be  calculated,  the  use  of  the  power  balance  as 
a  criticality  condition  would  result  in  an  expression  from  which  one 
could  obtain  the  entire  time -history  of  the  debond. 
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Figure  28.  Plot  of  strain  energy  versus  debond  area  for  end  loaded 
rocket  propellant  grain  debonding  axisymmetrically  from 
rigid  casing. 
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Figure  29.  Critical  pressure  versus  crack  length  for  bond  failure 

where  s  -  crack  length  and  b  =  radius  of  propellant  grain. 
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Whereas  the  application  of  the  pcwer  balance  in  predicting  the 
fracture  threshold  is  reasonably  straightforward,  the  calculations 
necessary  for  its  use  are  often  prohibitive  for  the  general  debond 
geometry.  As  indicated  in  the  preceding  paragraph,  this  is  especially 
true  when  the  problem  is  further  complicated  by  including  the  effects 
of  viscous  dissipation  as  one  should  when  considering  viscoelastic 
materials.  Thus,  as  in  the  case  of  the  spherical  flow  model  developed 
to  study  cohesive  fracture  in  viscoelastic  materials,  a  simple  model 
has  been  devised  which ,  it  is  felt,  leads  to  at  least  representative 
results  for  the  debonding  process  insofar  as  viscoelastic  dependence 
is  concerned.  Since  'he  analy." '  this  model  demonstrates  clearly 
many  of  the  principles  that  h»  a  discussed  in  this  section,  it 
is  appropriate  to  include  the  analysis  at  the  present  time, 
i  Consider  a  linear  viscoelastic  beam  of  infinite  length  which  is  bonded 
to  a  rigid  substrate  along  its  entire  length  with  the  exception  of  a 
portion  of  length  2a^  which  has  debonded.  A  stationary,  rectangular 
Cartesian  coordinate  system  x,  y,  z  is  oriented  as  shown  in  Figure  30  so 
that  the  origin  is  at  the  center  of  the  debond  and  the  x-axis  coincides 
with  the  axis  of  the  beam.  At  time  t  =  0,  a  pressure  q(t),  where 
q(t)  5  0  for  t  <  0,  is  applied  within  the  debond  causing  it  to  propagate 
symmetrically  at  some  time  t^  i  0  so  that  the  edges  of  the  debond  are 
at  |x|  =  a(t)  for  all  times.  Clearly  for  t  z  the  relationship 
a(t)  =  aQ  is  satisfied.  Neglecting  inertia  effects  and  using  beam 
theory,  it  is  a  simple  matter  to  solve  for  the  beam  deflection  w  and 
the  stress  6  and  strain  e  within  the  beam.  Because  of  the  viscoelastic 

X  X 

nature  of  the  beam  material,  these  quantities  are  functions  of  time  and 
are  given  by 


w(x,  t) 
6x(x,t) 
ex(x,  t) 


Ut; 

4V  J  L  J 


I3x2  -  a2 (t) ] 

y 


(82a) 
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where 


L'ht)  -  l_1[qCS)  SHcrp(S)] 


q(t)Dg  * 


n* 


3Dcrp(t'T) 

5T&-q(’)dl 


(83) 


and  I  is  the  moment  of  inertia  of  the  beam  cross-section  about  the  y-axis 
i  i 

and  L  denotes  the  inverse  Laplace  transform.  In  equation  (83),  the  bar 

over  a  function  denotes  the  Laplace  transform  of  that  function,  S  is  the 

transform  parameter,  Dcrp(t)  is  the  creep  compliance  of  the  material  and 

D  is  defined  as  D_=D_(0+). 
g  g  crpv 

As  a  criterion  for  debonding  of  the  beam,  an  application  is  made  of 
the  thermodynamic  power  balance  which,  neglecting  the  kinetic  energy, 
can  be  written  as 


2D 


SE 


(84) 


where, 

for  the 

present  problem, 

• 

I 

"  2af 

t>L 

di.  dx 

(85a) 

• 

d 

H2 

3.(t)  t 

C  f 

3e  (x,  x) 

F  +  2D 

"  2at 

T  dz 

/  *  J 

6  (X,  t)  - dr 

(85b) 

«^h/2 

0  -oc 

•  3x 

• 

a 

ait) 

• 

SE 

'  2af 

Jo 

■  2^aa(t) 

(85c) 

In  Equations  85  h  is  the  depth  of  the  beam,  the  dot  denotes  differen¬ 
tiation  with  respect  to  time  and  ya  is  considered  to  be  constant  with 
respect  to  time.  Also  it  should  be  noted  in  Equations  85a  and  85b  that  it 
is  necessary  to  integrate  with  respect  to  time  and  deband  area  or  beam 
volume  first  and  then  differentiate  with  respect  to  time  to  take  into 
account  the  change  in  debord  area  and  beam  volume  with  time.  Substi¬ 
tuting  from  Equation  82a  '.nto  Equation  85a  results  in  the  following 
expression  for  the  pwer  input  at  the  boundaries 
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i  ■  p  ar  r  *  L  11x2 ' a2(T,)2  L'1(,)l 

=  ijr  aMt  tw  >r  U»2w  -  ®2wi2  t'1wf  d» 
y  -oo 

(86) 

*  T5T;  ’(t)  It  u*2  •  a2(t)l2  L'l(t)I  * 

=  OT-  q(t)  a4W  a(t)  L_1(t)  +  jfp  q(t)  a5(t) 
yiy  y 

.  igi  j"  Bill  [a2(t)  -  a2(t)]2  L_1Ct)  dT  • 


Similarly,  substituting  Equations  (82b,  c)  into  equation (85b)  results  in 

a(t)  t  . 

F  +  2D  =  ^  f  d x  /  q(i)[3x2  -  aZ(T)]  ^  {[3*2  -  a  W]L'  (t)[  dx 


i-a(t)  f  q(t) [3a2 (t)  -  a2CO]  fj-  U3a2(t)  -  a2 (x]] L*1  (t)}  di 

V  J-oo 

(87) 

q(t)  [3x2  -  a2  (t)  J  {[3x2  -  a2  Cti  ]  L’ 1  (t] }  dx 

Jo 


'™~V 


gp  q(t)  a4(t)  a(t)  L  1(t)  +  j|p  q(t)  a^t) 


at 


.  1  ;(t)  Y  y  jq(T)[3a2(t)  -  a2(x)]j  [2  2(t)  -  a2^)]!"1^) 

Av  00 


di 


It  is  now  a  simple  matter  to  substitute  the  appropriate  expressions  into 
the  power  balance  equation  (84)  to  obtain 
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a(t) 


j  [3a2(t)  -  a2 (t) ] i.” 1  Ct)  ^  jq(T)[3a2(t)  -  a2(T)]  j  dr 

\ 

j*  ^aiil  (a2(t) .  »l(t)]2r1(T)  d,  -  72iy  ya 


(88) 


where  L  *(t)  is  given  by  equation  (83).  Equation  (88)  is  obviously 
satisfied  if  the  crack  does  not  run,  i.e.,  if  a(t)  *  0.  However,  setting 
the  bracketed  term  equal  to  zero  results  in  an  integral  relation  for  a(t) 
which  yields  yet  another  solution.  This  solution  is  the  time  history  of 
the  debond  growth. 

The  time  required  to  initiate  debonding  of  the  beam  for  a  particular 
loading  q(t)  cla  be  obtained  from  equation  (88)  by  noting  that  at  the 
time  of  initiation  a(t)  *  a^.  Setting  the  bracketed  term  in  equation  (88) 
equal  to  zero  and  a(t)  and  a(-r)  in  the  integrals  equal  to  a.g,  results 
in  the  expression 


91  v 
y’a 

a_ 


(89) 


where  again  L_1(t)  is  given  by  equation  (83).  Solving  equation  (89)  for 
t  yields  the  time  to  initiation  of  debonding  as  a  function  of  the  load  q(t) 
the  initial  debond  length  ag,  the  material  properties  and  beam  geometry, 
and  the  energy  necessary  to  create  new  adhesive  fracture  surface  y_. 

It  is  of  interest  to  evaluate  equations  (88)  and  (89)  for  two 
particular  loadings,  namely  a  pressure  applied  as  a  step- function  in  time 
and  a  pressure  applied  linearly  with  time. 


Step-Applied  Pressure*  3(f)  =  dgH(t); 

The  substitution  of  q(t)  *  qgH(t)  into  equation  (83)  results  in 

L_1w  ■  (so) 
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Substitution  of  this  expression  into  equation  (89)  then  gives 


Dcrp(T)  fi(T)dT* 


91  v 
y  Ta 

~*~r 


(91) 


where  6 (t)  is  the  Dirac  delta  function.  Hie  sifting  property  of  the  delta 

function  and  the  fact  that  D  _ (t)  =  0  for  t  <  0  results  in  the  reduction 

crpv  ' 

of  equation  (91)  to 


181  y 
y  ’a 

- % - 


(92) 


for  all  t  >  0.  Thus,  one  deduces  that  debonding  is  initiated  instantly, 
that  is,  at  t  =  0,  if  q^  is  greater  than  or  equal  to  the  qQ  defined  by 
equation  (92).  If  this  is  the  case,  then  the  time-history  of  the  debond 
for  t  >  0  can  be  found  from  equation  (88) .  Substituting  the  appropriate 
expressions  into  this  equation  and  setting  the  bracketed  term  equal  to 
zero  yields,  after  a  series  of  manipulations,  an  integral  equation  from 
which  one  can  solve  for  a(t)  for  t  >  0,  This  integral  equation  is  given 
by 


Hgj  2[3a2(t)  -  aQ2]2  -  3[a2(t)  -  a02]2j 

-  8  /  D 

JO  ' 


(93) 


.  1441  y 

(T)[3a2(t)  -  a2(x)]a(T)a(T)dT  = 


Since  what  is  of  primary  interest  is  the  initial  velocity  of  propagation  of 
the  debond,  equation  (93)  is  differentiated  twice  with  respect  to  time  to 
arrive  at 


a(0 

afiT 


2  n  f") 
crp 

15D  -  W'"Tt) 
g  crpv 


(94) 


where  it  has  been  assumed  that  a(t)  and  a(t)  are  both  nonzero  quantities 
for  t  >  0.  Thus,  since  debonding  is  initiated  at  t  -  0,  the 
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initial  velocity  of  propagation  is  given  by 

.  zy*> 

“o  05) 

Thus,  it  is  seen  that  the  simple  beam  model  proposed  does  not  predict  a 
smooth  transition  during  the  acceleration  from  zero  debond  velocity  at 
t  *  0’  to  the  initial  debond  velocity  at  t  ■  0+.  Also  it  is  noted  that 
Equation  94  is  not  entirely  satisfactory  since  for  larger  times  the 
velocity  of  propagation  of  the  debond  could  become  negative  for  certain 
materials. 

Constant  Pressure  Rate,  q(t)  «  qptH(t) 

Substituting  q(t)  ■  q0tH(t)  into  Equation  .83  leads  to  the 
result 


■  %Dcrp(1)M 


(96) 


where 


D,rp(1)M  *  jf  Dcrp(T)dT 


(97) 


Using  Equation  96  with  Equation  89,  the  time  to  the  initiation  of 
debanding  can  be  found  by  solving  the  equation 


D  ^  (O 
crp  v  - 


9I„  Y, 
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q0  a0 


(98) 


where 


D  (2) 
crp 


(t) 


Jo  Jo 


(c)dc  dT 


(99) 


Letting  the  solution  of  Equation  98  be  given  by  t  «  tf,  the  time-history 
of  the  debond  for  t  >  tf  can  be  found  from  Equation  88.  Substitution 

of  q(t)  *  qgtH(t)  and  Equation  96  into  Equation  88  and  subsequent 
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manipulations  loads  to  the  following  integral  equation  for  n(t): 
t. 


ft 


J?[3a2(t)  -  a2(r)]2  -  3[a2(t)  -  a  (t)]  |  Dcrp 


t  rn 

4  j  (3a2(t)  -  a2(,)]Ta(,)a(t)  DcrpUJ(t)  dt  - 


(100) 


in  arriving  at  ^nation  100  which  is  true  for  t  >  t.  the  bradcets  in 
Equation  (88)  were  set  equal  to  tero  since  for  t  >  tf,  a(t)  +  0  Differen 
tiating  Equation  100  with  respect  to  time  and  dividing  through  by 
a(t)  results  in  the  equation 

2D  ^(t)  a2(t)[a(t)  -  ta(t)] 
crp 


+  <i(t)  e  ja(3a2(t)  -  a2(r)]  -  3[a2(t)  -  a  (^j^crp^^ 

*o  l 

-  6a(t)  J  tD^^CO  a(T)a(t)dt  =  ° 


(101) 


which  upon  taking  the  limit  t  -►  t£  yields 


a(tf  ) 


D 


CD 


-  6Dcrp^J(tf) 


crp 


(tf) 


uF 


(102) 


.gain,  it  is  seen  that  the  simple  beam  model  does  not  predict  a  smooth 

transition  during  the  acceleration  fr«  zero  debond  velocity  at  t=£ 

,  ,  -a.  „*  «.  _  *  Also  i£  times  to  raiiure  aie 

the  initial  debond  velocity  at  t  -  tf  .  Also, 

Short,  Equation  102  may  predict  negative  initial  debond  velocities  where 
it  is  clear  that  these  values  must  be  positive.  For  example,  ta  ing  e 
limiting  case  of  a  glassy  material,  the  time  to  the  initiation  of  debon  ing 

is  given  by 


1RIy  Ya 

XT- 


(103) 
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t 


while  the  initial  velocity  of  propagation  of  the  debond  becomes 


(104) 


The  rather  unsatis factory  nature  of  the  results  obtained  for  the  model 
of  a  simple  beam  in  which  inertia  effects  have  been  neglected,  has  led 
to  the  consideration  of  a  model  which  includes  the  rate  of  change  of  the 
kinetic  energy  in  the  power  balance.  The  analysis  of  this  model  is  now 
proceeding.  It  is  thought  this  will  lead  to  more  satisfactory  results 
as  it  has  in  the  case  of  the  spherical  flaw  model. 
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V.  ESTIMATES  OF  TIME  TO  FAILURE  AND  INITIAL 
PROPAGATION  VELOCITIES  IN  BURNING  FLAWS 

In  order  to  obtain  first  order  estijnates  of  the  time  to  failure  and 
initial  crack  velocity  for  a  burning  flaw  in  a  propellant  grain,  it  has 
been  necessary  to  depend  on  existing  models  for  spherical  flaws  with 
slight  modifications.  For  the  case  of  the  time  to  failure,  it  was  assumed 
that  the  increase  in  radius  of  a  spherical  flaw  due  to  burning  was  linear 
in  time.  This  essentially  limited  the  analysis  to  problems  in  which  the 
average  pressure  in  the  spherical  flaw  remains  constant.  A  constant 
pressure,  as  discussed  in  Section  III,  is  not  realistic  unless  the  time 
to  failure  is  quite  small  since  the  burning  out  of  a  flaw  decreases  the 
pressure  an  order  of  magnitude  in  less  than  500  milliseconds  for  typical 
flaw  geometries  and  propellants.  The  analysis  to  determine  initial  pro¬ 
pagation  velocities  are  also  subject  to  the  assumption  of  constant  pressure. 
Discussion  of  the  time  to  failure  and  initial  velocity  are  given  in 
the  following  sections. 

FAILURE  TIME  FOR  BURNING  SPHERICAL  FLAWS  IN  PROPEM-ANT  GRAINS 

The  criticality  of  a  burning  spherical  flaw  is  considered.  A  solu¬ 
tion  is  obtained  by  use  of  Griffith's  crack  criterion  modified  for 
a  linearly  increasing  radius.  The  increase  in  radius  is  c^sumed  to  be 
due  to  burning  of  the  material.  The  analysis  considers  the  time  to  failure 
assuming  that  the  pressure  in  the  flaw  remains  constant.  This  analysis 
results  from  a  modification  of  Griffith’s  crack  criteria  by  Williams 
to  include  the  effects  of  different  geometries  and  viscoelasticity. 

Stated  mathematically  this  modified  criteria  is 

.  .  1(2  Ercl« 

-  p - 

c 
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where 

a 

K 


£rel«  ' 
YcCt)  - 


size  of  critical  flaw,  i.e.,  radius  for  spherical  flaw 

geometry  factor,  i.e.,  K  *  1.27  for  spherical  flaw 

relaxation  modulus ,  which  for  a  viscoelastic  material 
is  a  fuictic  1  of  time. 

surface  compliance  energy,  which  is  also  a  function 
of  time. 

internal  pressure 


The  variation  of  both  Erg^  (t)  and  rc(t)  with  time  can  be  approximated 
by  relations  of  the  form 


;rel«  * 

BtY 

(106) 

ycM  * 

6t* 

(107) 

where  B,  <5,  y,  and  4>  are  constants  to  be  determined  empirically. 

If  erosive  burning  is  neglected,  the  combustion  rate,  r,  of  a 
solid  propellant  may  be  assumed  to  follow  the  empirical  relation 

r  =  CPca  (108) 

where  r  is  the  recession  rate  of  the  burning  wall  and  C  and  a  are  cons ta.it s 
to  be  determined  empirically  for  a  particular  fuel.  From  Equation  108 
it  can  be  seen  that  the  assumption  of  constant  burning  rate  is  equivalent 
to  assuming  the  pressure,  Pc,  is  also  constant. 

From  Equation  108  one  can  also  determine  the  radius,  a,  of  a  flaw 
at  time  t,  assuming  linear  burning,  as 

a  =  aQ  +  CPcat  tlC9) 

ihere  a q  is  the  initial  flaw  size  at  time  t  =  0. 

The  flaw  size,  a,  can  now  be  eliminated  between  Equations  105  and  108 
to  obtain  the  final  result 
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(110) 


Equation  110  may  be  solved  interatively  for  t  by  use  of  Equations  106 
and  107  at  a  given  pressure  Pc  and  initial  flaw  size  aQ.  Figure  31  is 
a  plot  of  time  to  failure  versus  initial  flaw  size  at  ■various  pressures 
for  TP-H1011  propellant. 

As  can  be  seen  from  the  attached  figure,  for  all  pressures  investi¬ 
gated,  there  occurs  rapid  propagation  to  fracture 


INITIAL  PROPAGATION  VELOCITIES  FOR  BURNING  SPHERICAL  FLAWS  IN  PROPELLANT 
GRAINS 


In  the  preceding  section,  we  considered  the  time  to  the  initiation 
of  fracture  of  a  burning  spherical  flaw.  This  time,  tQ,  was  plotted 
versus  the  initial  flow  radius,  a  ,  for  several  different  pressures  in 
the  flaw.  Again  using  Williams1  $,19)  ^el  ^  a  spherical  flaw,  the  frac¬ 
ture  propagation  velocity  at  the  instant  fracture  takes  place  can  be  given 
as  a  function  of  the  time  to  failure.  This  result  is 


k  (2)  rt-j 

~rel  ltJ 


ss: 


Erel  (0  d5  dt 


(nib) 


Thus,  picking  the  time  to  failure,  tQ,  for  a  particular  initial  flaw  radius, 
a0»  from  the  graph  in  Figure  31,  this  value  of  t  may  be  used  in  Equation 
111  to  obtain  the  initial  fracture  velocity  a(tQ).  The  results  are  shown 
in  Figure  32,  where  a  graph  of  the  initial  fracture  velocity  versus  the 
initial  flaw  size  is  plotted  for  two  different  internal  pressures.  In 


Figure  32  are  also  shown  the  initial  flaw  sizes  for  the  two  pressures  at 
which  the  initial  fracture  velocity  exceeds  the  burning  rates.  For  initial 
flaw  radii  greater  than  these  values,  fracture  is  predicted. 
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CRACKING  RATE 
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VI.  CONCLUSIONS 


The  results  of  this  study  have  indicated  that  extremely  high  pres¬ 
sures  will  occur  in  burning  flaws  in  solid  propellant  grains .  The  smaller 
the  ratio  of  exit  flow  area  to  burning  surface  area,  the  larger  the  pres¬ 
sure.  Thus,  cracks  and  debonds  not  detectable  by  casual  visual  observation 
are  raosi:  likely  to  propagate  and  lead  to  rocket  motor  failure.  The  flaws 
detectable  by  visual  observation  will  be  of  equal  concern  only  if  they 
are  deep.  Debonds  such  as  those  occurring  in  the  Polaris  A-3  head- 
end  offer  same  experimental  verification  of  these  conclusions.  Debonds 
easily  detectable  by  ultrasonics  did  not  cause  motor  failure;  however,  on 
same  motors  without  detectable  flaws  failure  occurred.  The  experimental 
data  showed  both  case  bum  through  and  increased  specific  impulse  which 
could  be  correlated  with  the  additional  burning  surface  and  increased 
burning  rats  that  would  occur  in  a  burning  debond. 

Typically  the  critical  time  for  crack  propagation  is  500  milli¬ 
seconds.  By  this  time  the  flaw  lias  burned  out  sufficiently  that  the 
overpressure  has  dropped  by  an  order  of  magnitude.  Thus,  if  the  crack 
has  not  propagated  in  the  first  500  milliseconds,  it  is  doubtful  that 
it  ever  will. 

Limited  work  was  accomplished  to  determine  crack  propagation  ve¬ 
locities.  However,  utilization  of  a  spherical  flaw  model  provided  an 
estimate  of  the  initial  propagation  velocity.  This  estimate  indicated 
tha":  situations  could  clearly  develop  where  the  propagation  velocity 
exceeded  the  burning  rate  in  a  flaw. 


» 
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appendix  a 


numerical  solution 

Hie  simultaneous  solution  of  Equation  19  and  Equation  23  can  be 
approximated  numerically.  Equation  19  is  of  the  form 


P'  *  F(F,  n,  y) 

and  is  arr  table  to  approximation  by  the  simple  marching  method  which, 
assumes  that  the  derivative  is  constant  for  small  changes  of  the 
independent  variable.  With  this  approximation  one  can  write: 


Fi  + 


Ax  ¥' 


V 


(A-2) 


where  Ax  equals. the  step  size.  F  can  be  determined  at  each  step  by 
using  values  determined  at  the  previous  step.  The  only  variable  which 
is  unknown  is  n  and  it  can  be  approximated  by  numerical  integration  if 
F  is  considered  to  be  approximately  a  constant  equal  to  its  average  value 
over  the  short  interval  Ax  in  which  case  one  has  for  each  of  the  three  n 
integral  equations 
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computation  is  initiated  by  ccmputing  F£  from  Equation  A-l  with  f  determi 
from  Equation  25,  26,  or  27  for  the  geometry  under  consideration  with 


A-l 


^1  *  1*  nl  “  n0  "  ^V'1  ^  >"i  "  y0  Ca  function  of  the  given  geometry),  n 
is  then  computed  from  Equation  A- 3a,  A- 3b,  or  A- 3c  for  the  appropriate 
geometry  using  and  the  newly  computed  With  new  determined,  Fj  is 
confuted  ahd  so  on.  With  F^  and  known  at  each  step  and  y^  given,  other 
significant  parameters  of  the  flow  may  be  determined  such  as  the  Mach 
number  and  the  burning  rate. 

In  a  eiven  situation  the  material  properties,  flaw  geometry  and  chamber 

pressure  will  be  known.  From  these  known  parameters  an  initial  guess  at 

P  is  made  and  the  resulting  and  M  are  determined.  If  the  flow  is 
o  e  e 

choked,  P  is  increased  and  M  recalculated,  this  iteration  is  continued 
o  e  ' 

until  M  =  1.0.  If  the  flow  is  not  choked,  P  is  corrected  and  P  is  re- 
e  — —  '  o  e 

calculated;  this  iteration  continued  until  P  =  P  ,  . 

e  ch« 

The  difficult  portion  of  the  numerical  solution  lies  then  in  determining 
an  iterative  scheme  that  converges  rapidly  to  the  correct  PQ.  That  the  iter 
ations  converge  rapidly  is  paramount  since  each  guess  at  Pq  requires  that  at 
each  step  P- ,  n.  and  M.  be  recalculated  until  P  and/or  M  satisfy  the 

XXX  v  v 

required  boundary  condition.  In  a  typical  application  the  nwnber  of  steps 
will  vary  from  100  to  1,000  depending  upon  the  depth  of  the  flaw  and  the 
accuracy  required.  Appendix  B  contains  a  complete  listing  of  the  computer 
program  and  a  description  of  the  iterative  scheme  used. 

Because  the  solution  is  expressed  in  terms  of. the  pressure  at  small 
discrete  steps,  it  is  an  easy  matter  to  allow  for  change  in  geometry  due 
to  burning.  This  is  accomplished  as  follows: 

1.  A  pressure  distribution  is  determined  for  the  given 
initial  geometry  and  matched  to  the  boundary  conditions 
as  described  above. 

2.  A  small-time  interval  At  is  determined  such  that  the 
pressure  distribution  may  be  considered  approximately 
constant  over  At.  A  reasonable  At  is  the  time  that  it 
takes  yQ  to  increase  by  the  step  size  Ax. 

At  =  Ax/rbo  =  ax/CP” 


3. 


4. 


For  the  geometries  considered  the  new  y  dimension  as  a 
function  of  x  is  calculated  by 


rbi4t 


A  new  pressure  distribution  is  calculated  for  the  new 
geometry  and  matched  to  the  same  boundary  conditions 
as  before. 


(A-4) 


It  should  be  noted  that  as  burning  occurs,  the  assumption  that  y'«l 
is  no  longer  valid  in  the  region  of  the  tip,  since  this  region  assumes  the 
shape  shown  in  Figure  A-l  (if  no  cracking  occurs).  It  is  anticipated, 
however,  that  the  geometry  of  Figure  A-l  can  be  approximated  by  the  dotted 
line. 
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APPEMJIX  B 


A  complete  listing  of  the  computer  program  and  a  table  relating  the 
notation  used  in  the  analysis  to  that  used  in  the  program  is  located  at 
tiie  end  of  this  appendix.  Figure  1-B  illustrates  a  simplified  flow  chart 
of  the  program.  Notice  that  there  are  three  main  do-loops  labeled  1,  2, 
and  3.  Loop  3  serves  to  determine  the  umber  of  time  iterations  under¬ 
taken. 

Loops  1  and  2  both  determine  the  values  of  the  relevant  parameters 
(i^e.  P(I),  V(I),  etc.)  at  each  S(I)  for  a  given  PO.  Both  loops  contain 
a  series  of  IF  statements  which  direct  control  out  of  the  loops  to  a  state¬ 
ment  which  determines  a  new  PO  if: 

1.  P(I)  falls  below  PC  before  S(I)  =  L. 

2.  V(I)  rises  above  1.0  before  S(I)  =  L. 

3.  Z(I  +1)  becomes  negative. 

Both  loops  1  and  2  are  followed  by  a  series  of  IF  statements  which  deter¬ 
mine  if  a  correct  boundary  condition  has  been  met  (i.e.  P(N)  =  PC  and 
V(N)  £  1  or  P(N)  >  PC  and  V(N)  =  1)  and  if  so  directs  control  to  a  PRINT 
statement  and  if  not  to  a  statement  which  determines  a  new  PO.  Both 
loops  are  preceded  by  an  IF  statement  which  determines  if  the  maximum  num¬ 
ber  of  allowable  iterations  has  been  achieved  and  if  so  directs  control 
to  a  PRINT  statement  which  declares  that  the  iteration  is  not  converging 
rapidly  enough. 

The  object  of  having  two  DO  loops  which  are  so  similar  is  that  loop  1, 
in  general,  differentiates  the  choked  from  the  non-choked  flows,  or  in 
other  words,  loop  1  attempts  to  match  the  boundary  condition  that  P(N)  =  PC 
and  V(N)  £  1,  If  the  flow  is  either  choked  or  marginal  (i.e.  "just  about" 
choked)  then  loop  1  gives  a  "rough"  approximation  to  PO  and  control  then 
moves  down  to  loop  2  for  final  adjustment. 

Both  the  loop  1  and  loop  2  PO  correction  statements  are  based  on  the 
assumption  of  direct  proportioning,  or  in  other  wards,  it  is  assumed  that 
if  a  change  in  PO  of  AP0.  produces  a  change  in  P(N)  or  V(N)  of  AF(N),  then 

apo2  AP(N)2  AV(N)2 

A^  =  ATflD"  =  aV{NJY  (b_1) 


where  AP(N)^  or  AV(N)2  can  be  calculated  by  the  boundary  condition  one 
is  attempting  to  match.  In  reality  direct  proportioning  of  this  kind  is 
not  observed  and  one  must  reduce  the  amount  of  correction  aPC>2  by  dividing 
it  by  an  abritravy  constant  which  is  determined  by  the  boundary  condition 
being  met  and  the  size  of  the  correction.  For  instance  exit  Mach  number, 
V(N),  is  more  sensitive  to  PO  change  than  exit  pressure,  P(N),  and 
therefore  if  one  is  attempting  to  match  Mach  lumber  the  direct  proportion 
correction  PO  may  be  divided  by  10.0,  whereas  if  pressure  is  the  attempted 
match  the  division  constant  may  be  as  small  as  2.0  or  3.0. 

The  necessary  read- in  for  the  program  involves  three  sets  of  parameters; 
propellant,  geometry,  and  computer.  The  propellant  and  computer  parameters 
are  listed  at  the  beginning  of  the  program.  The  geometry  parameters  are 
Y(I),  YP(I)  and  S(I).  They  may  be  read-in  by  either  assigning  a  value  of 
Y(I)  and  YP(I)  to  each  value  of  S(I)  or  by  defining  a  functional  relation- 
slip  between  Y(I),  YP(I)  and  S(I).  In  either  case  one  must  keep  in  mind 
that  the  program  was  developed  assuming  that  YP(I)  is  small. 
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Figure  1-B  Simplified  flow  chart  of  numerical 
analysis  computer  program 


TABLE  B-l 


Nomenclature  and  Units  Used  in  Computer  Program 


Analysis 

Computer 

Units 

Br 

BR(I) 

none 

C 

C 

£t^n+1/sec-lb: 

D 

D 

lbm/ft3 

f 

FC 

none 

G2 

F 

(sec-lbf/lbm)‘ 

H 

HBR 

none 

M 

v(i) 

■none 

n 

XP 

none 

P 

p(i) 

psi 

Po 

PO 

psi 

Pch 

PC 

psi 

P 

Z(I) 

none 

V' 

ZP 

1/in 

k 

R 

ft-lbf/lbm°R 

Re 

RE(I) 

none 

To 

TO 

°R 

At 

DT 

sec 

X 

SCD 

in 

AX 

H 

in 

Y 

YCD 

in 

Y 

YP(I) 

none 

6 

(34 

none 

C 

X(I) 

none 

y 

u 

lbf-sec/ft^ 

n 

T(i) 

none 

To 

TOO 

none 
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HAP  0017-0I/29-Q7S30 

STaRT»0Q607i  t  PROG  S1ZEU/0)»3005/I3uA 

DIMENSION  Yi2000)  |YP(2000I  ,S(2t*00)  ,Z(2000J  ,x(200o»  • 

•  Y<2000| 

DIMENSION  P ( 2000 1 • T ( 2000 ) ,  Re <2000) . d« ( 2000 ) , G ( 2UnO I 
read  n7,xp,rito,u,c*o»<?m 
'  H7  FOKhaTI  IOcU.S) 

C  PROPELLANT  PARAMETERS 

c  xP^OURMNG  rate  exponent, 

c  RaGAS  CONSTANT  FT  LBF/L3M  R 
C  TO“flURM I  Mg  TEMPERATURE  R 
C  U" V  1 SCOS 1 T Y  LBf  SlC/fT**2 

c  C*0UKM1HG  RfcTE  COEFFICIENT  FT**(2XP*1)/SEC  LBF**xP 
c  D»PROPELLaNT  DENSITY  Lflh/FT »#3 
C  6M*R AT  1 0  or  SPECIFIC  NEATS 

READ  N8,  D1FiERRCRiHiEiMAX»M|N»NA 
58  F0RMAT<3E10«5.SI I0> 
c  COMPUTER  PARAMETERS 

C  DIF»ALLOAA0LE  ERROR  BETWEEN  ACTUAL  Mach  no  and  MACH  1*0 
c  ERROR-ERROR  IN  PRESSURE  BETWEEN  CHANBtR  and  exit 
C  H»STEP  SIZE  IN  INCHES 

C  L«NO  OF  STEPS  PRINTED  IE  PRINT  EVE«Y  L  TH  STEP 
C  MAX«HAXIMUM  NO  OF  ITERATIONS  ALLOf  Ed 
C  M-NO  OF  TIME  INCREMENTS  TAKEN 
c  n»no  or  steps  Taken 

c  NA-NO  OF  STEPS  PRINTED  -  IE  NA  STE^S 
C  INITIALIZE  Tlr.E  STEP 
DT  ■  0. 

C  THIS  LOOP  DOES  THE  TIME  STEPS 
00  ]1  J  *  |,  H 

c  INITIALIZE  ITERATION  COUNTER 
N  1  «  1 

C  INITIAL  PO  g°ESS 
P0«PC*1 , 

PORS  -PC 
SR"0» 

GO  TU  2a 
28  CONTINUE 
opos»po-pors 
OS  cS(I»-SR 
IF(dS) 100, 1 U 1 ,100 
101  DS«H 
100  CONTINUE 
SRBS ( I ) 

PORS*PO 

IF(NI-6)S<lfS<t,S3 
55  SDIV-2.S 
GO  TO  57 
53  Sd! |  • 

57  continue 
C  CORRECTED  Po  GUESS 

PO-POMS»iDPOS/SOl V )  * (S (N  )  -  SR  */DS 
26  CONI  I NUE 
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T«0  ■  U*C*(P0*  1  HS*  J  •  XP-1  .  ) 

TCI)  ■  0*C«SP0*kR‘»»)  «♦  C  XP-l  •  )  . 

Z  €  1  )  *  1  • 

F«  I  2  •  +  I 6M* 1,1  •  TO  •  R»/(GH  •  32*2) 

c  upgrade  iteration  counter 

NI  *  MI* I 

C  CHECK  IF  MAX  ITERATIONS  HAVE  BEEN  REACHED 
IF  INI -MAX >58, 58,29 
58  CONTINUE 

C  THIS  LOOP  DETERMINES  VII)  AND  P  f I )  AT  EACH  STEF 
DO  3  '*  I  *  1*  IJ 
F(I)  «  PO  •  2!  I  ) 

’ F I  I “N ) SO ,22,22 

50  CONTINUE 

C  EXIT  LOOP  IF  PH)  FALLS  BELO«  Pc 

I F  (  P  (  I  )*ERR0R-PC)2*I,22»22 

22  CONTINUE 

XII  RSQnTIZI  |  l«ZI  I  )+F«TI’!)  *T(  I  )  ) 

V(I)=(X<I)-Z(I  ))/(/( 1  )• ( gh**  |  « )  ) 

It  IRNISI  ,23,23 

51  CONTINUE 

C  EXIT  LOOP  IF  V.I)  RISES  ABOVE  1,0 
IF (VI  I >-l »  123,29,29 

23  CONTINUE 

RE  I  1  >»IYl  1  /  »T  I  I  1 *P0*29 •  l/(U»32>2) 

C  CHECK  IF  FLO*  IS  TURBULENT  Of  LAMINAR 
IF  (REI  I)  -  2200* )  20,  20,  21 

20  FC  B  96,  /RE  II) 

GO  TO  27 

21  FC  «  *316/  Re(I)  ,25 
27  CONTINUE 

BRI  I J«THO*ZI I >*«XF*SQRT(RE< I ) ) / T (  I  ) 

52  CONTINUE 

IFIN-I ) 1O81IO8, 105 
105  CONI lNUE 

HBR«1  ,/12.«IBR  <D*|,  11*1.  /(  2,»5,*«  I  OR  (  |  )/<«•)  ) 

B2  ■  F  •  TCI)  /  XU) 

89  *  <GH*  XI  I  )  )  /  I  (XII)  -  gm  ♦  ZIP)  •  Vim 

c  use  this  m  card  for  a  dfbonb 

HlBlA(ll-ZII))*lPC*(I. *HHR ) /2,*YPI I ) ) 

C  USE  THIS  3 1  CARO  for  a  crack 

BI*IM  I  )-Tl  I  )  )«tFC*HBR*YPt  1  )  I 
C  USE  THIS  Bl  CaRO  FOR  A  CIRCULAR  Ut^OND 

Bl»l XI  I 1-ZI I ) » • l  PC*  I 1 •*H)/2,*yP( I 1-YI I l/IRO-St 1 1  )  I 
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c  use  this  b3  caru  for  a  oeqond 

B3«THO*Z ( I ) ••XP-YP I  I  )  «T ( I  ) 

C  USE  THIS  83  CARD  FOR  A  CRACK 

B3*2  •  •  THO«Z(  II»*XP-YPm*Tt !  J 

c  use  this  &3  card  for  a  circular  debond 

B3"Tho«Z<  I  )*«XP-YPI  I)iTm*?i;»»Hl)/|hO-S(i  I  I 

c  use  this  as  card  for  a  circular  o*$ond 

B5*Z»  I  )  t-Y  (  I  ) / I RO“S (  I  Ut(CH-t  •  )/6H 

c  use  this  zp  card  for  a  oebond 

ZP  •  BR  • CB I  ♦  B2  •  B3 )  ' 

c  use  this  zp  card  for  a  crack 

ZP  ■  BR  *(Bt  ♦  82  •  23) 

c  use  this  zp  card  for  a  circular  dibond 
ZP»<U1  ♦ab*B2*B3  S*B*< 

Z(l«l|  .  Zill  ♦  H*ZP 
JF<  Z ( I ♦! > Vz ( I ) )2R,2M»2S 
25  CONTINUE 

C  USE  THIS  8  CARD  FOR  A  DEBOMO 

B  «•  121)41)4  Z( J ) XP  4  * 

C  USE  THIS  R  CARD  FOR  A  CRACK 
Bs2t«  (  Z(  1  **2(  l4llU»XP 
C  USE  THIS  ft  CARD  FOR  a  CIRCULAR  OEBOND 

B» ( Z ( I ) *Z ( I  4  I ) )4»xP*  IRO-H4I4H/2*  ) *8 
C  USE  THIS  T  < I ♦ 1 »  CARO  FOR  A  DEBOND 

T 1 1  ♦  1  )»CThO/YU*II  )»IH*B/2«*4XP4Y0) 

C  USE  THIS  TII+1)  card  FOR  A  CRACK 

T(  I*  I  )»lTHO/Yl  !  +  m«(H*B/2t**AP4Y0) 

C  USE  THIS  T ( ! ♦! |  CARD  FOR  A  CIRCULAR  DEBOND 

?.<!♦!  I  «  THO/ <  Y  C 1  ♦  2  )*IR0-SI  1*1  )  )  ) • t H»e/ < 2 » • • Xp ) 4 Y <  J  I  *Rq  > 
108  CONTINUE 

VC  1 1  "SQRT  c  V  C  I.)  I 
36  CONTlNUC 

THE  FOLLOWING  «|  if  STATEMENTS  DETERMINE  IF  A  CORRECT 
BOUNDARY 

condition  has  been  satisfied 

IFCV  f  N  I  *0  l  F  •  >60,60,59 
80  IF<P(N)4’ERR0R-PC)V0iS6,9I 

91  !FlABS!P(N)-pC>-tRR0R)56,S6,?2 

92  IFC  ABS<  VCNI-J  •  > -D 1 F J b6 , 56 , 9 j 

C  THIS  if  STATEMENT  DETERMINES  WHICH  BOUNDARY  CONDITION  TO 

c  attempt 

C  TO  MATCH 

93  IFCAtiS(P!N)-PC)/PC-ABS(I #-V(N) ) ) 90, 90 ,88 
59  DPO-PO/IUO. 

60  TO  89 
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90  0P0»(PC-P(N) )/10* 

GO  TO  09 
80  UPO»  PO/lOO. 

8?  continue 
c  cokrectcd  po  guess 

PO»PO*DPO 

EX0P-0P0/10, 

NV*  I 
NP- 1 

55  CONTINUE 
VR*  V  (  I  > 

PR*P {  I  ) 

c  upgrade  iteration  COUNTER 

NI  a  N  I  ♦  1 

C  CHECK  IF  MAX  ITERATIONS  HAVE  BEEN  REACHED 
IFCNI-  MAX ) 66 , 86 , 29 
86  CONTINUE 

B  *  U. 

ThO  *  0*C«(Po« INN* ) *« (XP*I  •  ) 

Tin  a  o*c«(po#m» )  ••ixp-i.i 

l  I  1  I  a  I  . 

F“  I  2  •  •  C  <jM*  I  •  )  •  TO  •  RI/lGM  •  32.2) 

C  THIS  LOOP  OeT  eRm  I  Nj_S  VID  AND  Pjl)  AT  EACH  STEP 
DO  37  I  a  1  tN 
PCI )  «  pO  «  ZCI) 

I  F  C  I  **N  )  6  I  ,62,62 
6!  CONTINUE 

C  EXIT  LOOP  Ip  PH)  FALLS  BELOW  Pc 
IFCPC  I  )  ♦ERR0K-PC)2«C  ,62,62  1 

62  CONTINUE 

X ( I ) "SORT ( Z ( I )*ZC I l+F*T(  I  »*T<  1 1  I 
V{1)"CX(I)»Z(I))/(2C I )*(GH-J .  )  ) 

I  F  I  I-N>63,6H,6i» 

63  CONTINUE 

C  EXIT  LOOP  If  VCI)  RISES  ABOVE  1.0 
1FIVII)-I»)6*|,6II*2<I 
6S  CONTINUE 

R£C  I  IsfyC  ]  )*T(  I  )  •  P  0  *  2  R  »  )/(U*32*2) 

C  CHECK  If  FLOW  IS  TURBULENT  OF  LAMINAR 
IFlKtC I ) "220 J* ) 65 , 65 ,66 

65  FC J V6 , /RE  I  I ) 

60  TO  67 

66  FC»*3 1 6/Rf  I  I ) ••  .25 

67  CONTINUE 

BRC I  )aTnO.Zm#*XP»SQP  (RE  (  I  J  >/Tf  I  ) 
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|f <N- I ) 106 c 1 06 t |07 
107  CONTINUE 

»2  «  f  •  id)  /  mn 

BI  ■  ( GH*  x  I  I  1  »  /  (<X(I)  -  SM  •  Z(ln  •  7(1,) 

c  use  this,  st  card  for  a  debond 

B 1  •  (  X  (  I  )“Z(  I  )  I  •  1 F  C  •  1  l .♦HBR>/2.*YP(  I  )  ) 

c  use  this  bi  card  for  a  crack 

BI • ( X ( | J  —  Z 1 1 ) ) • ( FC*HBR* YP ( I  i  ) 

c  use  this  bi  card  tor  a  circular  debond 

0!«(X(1 I -z  « I ) )*(FC*( 1 .♦H)/2.*YP( I )-T ( I J / «  KO-S ( I )  }  ) 

C  USE  THIS  B3  CARD  FOR  A  DeBOND 
B3«TH0»Z< I > •• XP-VP I  I ) *T I  1 1 
c  use  this  <u  card  for  a  crack 

B3*2 • •  THO«Z{ I)#«XP-YP( I )*T(  I  ) 

C  USE  THIS  63  CARD  FOR  A  CIRCULAR  DE*OWD 

B3»TH0*Z<!l»#Xp-YP(I  )#T(|  |ty|  |  lull  )  /  (  RO-S  (  j  )  ) 

C  USE  THIS  as  CARD  FOR  A  CIRCULAR  DEBOND 
85*2  I  I  I *T( I l/IRO-SC I ) >«(GH-| • )/GM 

c  use  This  zp  card  for  a  debono 

ZP  »  BB  • C  8 1  ♦  82  *  B3) 

C  USE  THIS  ZP  CARD  FOR  A  CRACK 
ZP  ■  BM  •(81  ♦  82  •  B3 I 

C  USE  this  zp  CARO  FOR  a  circular  debond 

ZP«(0|*b5*B2*B3)*BH 
2(1*1  I  «  Z(I)  ♦  H*ZP 
I F ( Z ( I ♦ I ) *Z I  I ) )2H,69,69 
69  CONTINUE 

C  USE  THIS  0  CARD  FOR  A  DEBOND 

B  *  ( Z ( I  *  I ) *  Z ( I  I  )••  XP  ♦  B 

c  use  this  b  card  for  a  crack 

B*  2  ♦  •  (  Z  (  II*Z(  1*1  n*«XP  *B 

c  use  this  b  card  for  a  circular  de&ond 

8a(Z(l)*Z(l  +  l)  J**xP«(R0-H»I+H/2,  )  *B 
C  USE  THIS  T ( I  *  I  I  CARO  FOR  A  DEBOnD 

T(  I  ♦!  )"(THO/Y  (  I  +  m»(H«B/2***XP*Y0) 

C  USE  THIS  T ( I ♦ I )  CARO  FOR  A  CRACK 

T( I ♦ I )«(Th0/y ( l ♦ I ) ) *(H*B/2. •*XP*y0) 

C  USE  THIS  T ! I ♦ J )  CARO  FOR  A  CIRCULAR  DEBONO 

T( 1*1 )«THO/(Y( 1*1 ) • ( KO”S ( 1*1  I) )*(H*B/(2»**Xp)*V ( i )*RoI 
106  CONTINUE 

V ( I )>SQrT< V( I  I  ) 

37  CONTINUE 
DP»PIN»-PR 
DVbVIn)-VR 
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0RC«A8S(P(N»-PC)/FC-ABS( I.-VCNJ  ) 

C  THE  FOLLOWING  4  IF  STATEMENTS  OETERMINE  IF  A  CORRECT 

C  boundary 

condition  has  been  satisfied 

IF(P(N)+ERHOR-PCj8  0,8l  ,81 

81  IF< V ( NJ-D1F-1 •) 82,56,83 

82  IF«ABS(PlN>-PC)-ERROR)S*.RAtfl7 
87  !F(ABS<V(N)-J')-0!F)St.,S6,8<i 

c  this  if  statement  determines  which  boundary  condition  to 
c  ATTEHPT 

c  to  match 

84  IF(AbS^P(N)-PC)/P(N|-ABS( 1 .-V(N» ) 180,80,83 

80  CONTINUE 
NP-NP+I 

CNP*A8S(0p/EXDP  1 
IF<nP-2)9S*95.94 
9S  PD  I  V  »3 • 

60  TO  99 
94  PDIV-CNP 
99  COHllNUE 

POR»PO-( (P(N)-PC»/POIV)*OPO/DP • 

EX0P»(P(N)-PC)/P0IV 

GO  TO  6S 

83  CONTINUE 
NV"N V* 1 

lF(NV-2)97, 97.96 

98  IF« A^SI V (N)-l . )-«3) I  04 , 103 , I  03 
|04  VO  I V*2  » 

GO  TO  98 

97  IF< ABSI V  «  N I - 1 . I -.3)  1 02, 103, 103 

102  V0lV«3. 

GO  TO  98 

103  V0IV*10*0 

98  CONTINUE 

POR»PO-( ( VINJ-l , )/V0IV)«0PO/DV 

85  CONTINUE 
OPO*POR-PO 
PO-POR 

GO  TO  55 
56  CONTINUE 
’ DR»D«*OT 

C  PRINT  DESIRED  PARAMETERS  WHICH  JOeNTIFY  INITIAL  GEOMeTRY, 
C  TIME  , 

C  PROPELLANT,  AND  PC 

C  PRINT  OESJREO  VARIABLES  ALONG  LENGTH  OF  FLAW  IE  V(I), 
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c  pllh  *<!»#  S(!) 

29  CONTINUE 
ORaQK+OT 

C  PRINT  DCS | RED  PARAMETERS  WHICH  IDENTIFY  INITIAL  GEOHeTKY, 
C  TIHE  , 

c  PROPELLANT#  ano  pc 
PRINT  72 

72  FORMAT!*  *i«lTERATION  unstable*! 

26  CONTINUE 

DTaH/IC»(POM9*f«)**XPM2»  ) 

C  COMPUTE  NEW  GEOMETRY 
DO  3b  I  a  1  f N 
Y(II*r(!!*H«2<n»»XP 
Y II > ■ Y i I »+H*Z( I»«aXP*2. 

35  CONTINUE 

00  36  I  a  ||N 

ypiii  a  mi«n-Y(m/  h 

36  CONTINUE 
39  CONTINUE 

STOP 
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APPENDIX  C 


Consider  the  plane  problem  of  a  Griffith  crack  of  length  2a  in 
a  thin  elastic  sheet.  If  a  uniform  pressure  acts  within  the 

crack,  the  symmetric  stress  distribution  in  the  vicinity  of  the  crack 
tip  is  given  by 


Ki 

4/2  7 rr 
Ki 

4/2t rr 

Ki 

4/2’ir 


(C-la) 
(C-lb) 
iC- lc) 


where  r  and  9  are  polar  coordinates  with  origin  at  the  crack  tip  as 
shown  in  Figure  C-l,  and  K,  is  the  node  I  stress  intensity  factor 
wh^ch  for  this  problem  can  be  written  as 


K,  -  PQ  /S 

Applying  the  energy  balance  as  a  fracture  criterion  leads  to  the 
criticality  condition 


cr 


(C-3) 


This  condition  can  be  put  into  the  more  familiar  but  equivalent  form 


where  y  is  the  cohesive  fracture  energy,  by  recognizing  that  the 
potential  energy  release  rate  with  crack  extension  is  directly  related 
to  the  stress  intersity  factor  Kj  ^17^. 

If  the  pressure  inside  the  crack  were  now  to  be  distributed  non- 
uniformly  along  the  length  of  the  crack,  it  becomes  important  to  assess 
the  effect  of  this  non-uniformity  upon  crack  instability.  A  convenient 
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way  of  doing  so  is  to  caapare  the  fracture  threshold  for  the  non¬ 
unifora  pressure  distribution  to  that  which  would  result  if  its  total 
(integrated)  pressure  were  distributed  along  the  crack  length  uni¬ 
formly.  This  comparison  entails  solving  for  the  stress  intensity 
factor  for  the  non-uniform  pressure  which  will  be  accomplished  by 
using  the  Green's  function  obtained  by  the  method  of  14iskhelishvili 

If  two  oppositely  directed  concentrated  forces  each  of  magnitude  P 
act  within  and  normal  to  a  crack  at  some  distance  b  from  its  center,  as 
shown  in  Figure  C-l,  the  stress  intensity  factor  is  given  by 


Recognizing  that  Equation  (C-5)  is  a  Green's  function  for  the  problem 
of  a  non-unifora  pressure  p(x)  applied  within  the  crack,  the  stress 

intensity  factor  for  this  problem  is  easily  obtained  in  the  integral 
fora 


Approximating  the  pressure  within  the  crack  by 


P(0  -  P0[l  -A(l-£2)m) 


(C-7) 


where  ?  =  x/a  and  X  and  m  are  parameters  chosen  to  K*st  fit  the 
actual  pressure  distribution,  one  obtains  the  stress  intensity  factor 
for  this  distribution  from  Equation  (C-6)  as 

Ki  ■  S*  t [i  • (*$■ « 

■  p0tra)!i  [  1  -  ?  B(”.  *  7>  |)]  (C-8J 
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In  Equation  (C-8),  B(a,B)  is  the  complete  Beta  function.  Upon 
expressing  Equation  (C-8)  in  tenns  of  the  average  of  the  applied 
pressure  p,  one  finally  obtains 


Kj  =  p(tra) 


1  -  (A/it)  B(m  +  £) 
1  -  (A/2)  B(m  +  1,  |) 


(C-9) 


where 

P  -  \  P(5)d?  -  p0  £l  -  j  B(m+l,|)J  (C-10) 

Note  in  Equation  (C-9)  that  if  A  =  0,  the  uniform  pressure  result  is 
recovered.  Hence  the  term  in  brackets  in  Equation  (C-9),  which  is 
graphed  in  Figure  C-2,  expresses  the  deviation  in  the  stress  intensity 
factor  due  to  deviations  of  the  assumed  pressure  (C-7)  from  uniformity. 
This  deviation  is  expected  to  be  typical  of  non-uniform  pressure  distri¬ 
butions  in  cracks  of  rather  arbitrary  geometry. 
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